[BioC] problem normalizeMedians (limma)

[Marcus Davy]

Hi,
Doe this include modifying normalizeMedianDeviations aswell for consistency.
I suspect normalizeMedianDeviations doesn't reach R's limit as quickly.

Marcus


MA <- new("RGList")

n     <- 200
genes <- 1000

MA$M <- matrix(runif(n*genes, -8, 8), nr=genes, nc=n)
MA$A <- matrix(runif(n*genes, 1,16), nr=genes, nc=n)

summary(c(MA$M))
summary(c(MA$A))

any(is.infinite(normalizeMedians(MA$A)))

normalizeMedians2 <- function(x)
{
  narrays <- NCOL(x)
  if (narrays == 1)
    return(x)
  a.med <- log(apply(x, 2, median, na.rm = TRUE))
  a.med <- exp(a.med - mean(a.med))
  t(t(x)/a.med)
}

# Subset
normalizeMedians(MA$A[1:5,1:5])
normalizeMedians2(MA$A[1:5,1:5])

# Entire dataset
any(is.infinite(normalizeMedians(MA$A)))
any(is.infinite(normalizeMedians2(MA$A)))

normalizeMedianDeviations2 <- function(x)
{
  narrays <- NCOL(x)
  if (narrays == 1)
    return(x)
  medabs <- function(x) median(abs(as.numeric(x[!is.na(x)])))
  xmat.mav <- log(apply(x, 2, medabs))
  denom <- mean(xmat.mav)
  si <- exp(xmat.mav - denom)
  t(t(x)/si)
}

normalizeMedianDeviations(MA$M[1:5, 1:5])
normalizeMedianDeviations2(MA$M[1:5, 1:5])

# Entire dataset
any(is.infinite(normalizeMedianDeviations(MA$M)))
any(is.infinite(normalizeMedianDeviations2(MA$M)))



On 10/19/05 1:19 PM, "Gordon Smyth" <smyth at wehi.edu.au> wrote:

> 
>>> Date: Mon, 17 Oct 2005 18:24:47 -0400
>>> From: Francois Pepin <fpepin at cs.mcgill.ca>
>>> Subject: [BioC] problem normalizeMedians (limma)
>>> To: bioconductor at stat.math.ethz.ch
>>> 
>>> Hi everyone,
>>> 
>>> I'm not sure what should be done about it, but I've been surprised with
>>> the following behavior with normalizeBetweenArrays:
>>> 
>>> (dat is an MAList here)
>>>> dat$A[1,1]
>>> [1] 8.796361
>>>> dat <- normalizeBetweenArrays(dat)
>>>> dat$A[1,1]
>>> [1] Inf
>>> 
>>> This is due to the fact that I've got almost 400 arrays here. The
>>> normalizeMedians(x) (which is called by normalizeBetweenArrays) has the
>>> following behavior:
>>> 
>>>     a.med <- apply(x, 2, median, na.rm = TRUE)
>>>     a.med <- a.med/(prod(a.med))^(1/narrays)
>>> 
>>> multiplying 400 relatively small numbers can easily reach R's limit as
>>> far as doubles are concerned.
>>> 
>>>> 6^400
>>> [1] Inf
>>> 
>>> Wouldn't it be better to be working in log space instead?
> 
> Yes, you're right. In limma 2.3.2 I will change
> 
>     (prod(a.med))^(1/narrays)
> 
> to
> 
>     exp(mean(log(a.med)))
> 
> which should make floating underflow very unlikely.
> 
> BTW, the "scale" method has been the default for normalizeBetweenArrays()
> for historic reasons, inherited from the older sma package. I'm going to
> change the default to "Aquantile".
> 
> Gordon
> 
>>> Francois
>>> 
>>>> sessionInfo()
>>> R version 2.1.0, 2005-04-18, x86_64-unknown-linux-gnu
>>> 
>>> attached base packages:
>>> [1] "methods"   "stats"     "graphics"  "grDevices" "utils"
>>> "datasets"
>>> [7] "base"
>>> 
>>> other attached packages:
>>>   limma
>>> "1.9.6"
> 
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/bioconductor


______________________________________________________

The contents of this e-mail are privileged and/or confidenti...{{dropped}}