[BioC] limma normexp background correction bug?

[James W. MacDonald]

Tobias Straub wrote:
> hi
> 
> i was playing around with various background correction methods using 
> limma and discovered that whatever normexp.method i call for the 
> "normexp" method, the result is exactly the same. in fact the more i 
> digged into it the more i believe that normexp is not executed at all.
> 
> looking at the code for backgroundCorrect it more looks like some junk 
> of code has been lost. i also checked the latest dev version of limma 
> (3.0.3), looks the same. instead of doing normexp the execution jumps 
> into some "rma" kind of thing that just leaves me confused.

You are getting confused because you miss the point that

switch(x,
foo = do.foo(),
bar = do.bar(),
foobar = , barfoo = do.barfoo(),
)

means that if the argument is foobar _or_ barfoo, you do the same thing.

So normexp *is* used, but it appears that if you pass in either method = 
"rma" or method = "normexp", you will always end up using method = 
"saddle" for normexp.fit().

Whether or not this is intentional or a bug is something only Gordon can 
decide.

Best,

Jim


> 
> any clues someone? any correct code around somewhere?
> 
> best
> T.
> 
> 
>  > backgroundCorrect
> function (RG, method = "subtract", offset = 0, printer = RG$printer,
>     normexp.method = "saddle", verbose = TRUE)
> {
>     if (!is.list(RG) && is.vector(RG))
>         RG <- as.matrix(RG)
>     if (is.matrix(RG)) {
>         method <- match.arg(method, c("none", "normexp", "saddle",
>             "neldermean", "bfgs", "rma", "mcgee"))
>         if (method == "normexp")
>             method = "saddle"
>         if (method != "none") {
>             for (j in 1:ncol(RG)) {
>                 x <- RG[, j]
>                 out <- normexp.fit(x, method = normexp.method)
>                 RG[, j] <- normexp.signal(out$par, x)
>                 if (verbose)
>                   cat("Corrected array", j, "\n")
>             }
>         }
>         if (offset)
>             RG <- RG + offset
>         return(RG)
>     }
>     if (is.null(RG$Rb) != is.null(RG$Gb))
>         stop("Background values exist for one channel but not the other")
>     method <- match.arg(method, c("none", "subtract", "half",
>         "minimum", "movingmin", "edwards", "normexp", "rma"))
>     if (method == "rma") {
>         method <- "normexp"
>         normexp.method <- "rma"
>     }
>     normexp.method <- match.arg(normexp.method, c("mle", "saddle",
>         "rma", "rma75", "mcgee"))
>     if (is.null(RG$Rb) && is.null(RG$Gb))
>         method <- "none"
>     switch(method, subtract = {
>         RG$R <- RG$R - RG$Rb
>         RG$G <- RG$G - RG$Gb
>     }, half = {
>         RG$R <- pmax(RG$R - RG$Rb, 0.5)
>         RG$G <- pmax(RG$G - RG$Gb, 0.5)
>     }, minimum = {
>         RG$R <- as.matrix(RG$R - RG$Rb)
>         RG$G <- as.matrix(RG$G - RG$Gb)
>         for (slide in 1:ncol(RG$R)) {
>             i <- RG$R[, slide] < 1e-18
>             if (any(i, na.rm = TRUE)) {
>                 m <- min(RG$R[!i, slide], na.rm = TRUE)
>                 RG$R[i, slide] <- m/2
>             }
>             i <- RG$G[, slide] < 1e-18
>             if (any(i, na.rm = TRUE)) {
>                 m <- min(RG$G[!i, slide], na.rm = TRUE)
>                 RG$G[i, slide] <- m/2
>             }
>         }
>     }, movingmin = {
>         RG$R <- RG$R - ma3x3.spottedarray(RG$Rb, printer = printer,
>             FUN = min, na.rm = TRUE)
>         RG$G <- RG$G - ma3x3.spottedarray(RG$Gb, printer = printer,
>             FUN = min, na.rm = TRUE)
>     }, edwards = {
>         one <- matrix(1, NROW(RG$R), 1)
>         delta.vec <- function(d, f = 0.1) {
>             quantile(d, mean(d < 1e-16, na.rm = TRUE) * (1 +
>                 f), na.rm = TRUE)
>         }
>         sub <- as.matrix(RG$R - RG$Rb)
>         delta <- one %*% apply(sub, 2, delta.vec)
>         RG$R <- ifelse(sub < delta, delta * exp(1 - (RG$Rb +
>             delta)/RG$R), sub)
>         sub <- as.matrix(RG$G - RG$Gb)
>         delta <- one %*% apply(sub, 2, delta.vec)
>         RG$G <- ifelse(sub < delta, delta * exp(1 - (RG$Gb +
>             delta)/RG$G), sub)
>     }, normexp = , rma = {
>         if (verbose)
>             cat("Green channel\n")
>         RG$G <- backgroundCorrect(RG$G - RG$Gb, method = method,
>             verbose = verbose)
>         if (verbose)
>             cat("Red channel\n")
>         RG$R <- backgroundCorrect(RG$R - RG$Rb, method = method,
>             verbose = verbose)
>     })
>     RG$Rb <- NULL
>     RG$Gb <- NULL
>     if (offset) {
>         RG$R <- RG$R + offset
>         RG$G <- RG$G + offset
>     }
>     new("RGList", unclass(RG))
> }
> <environment: namespace:limma>
> 
>  > sessionInfo()
> R version 2.9.2 RC (2009-08-17 r49303)
> x86_64-apple-darwin9.7.0
> 
> locale:
> en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8
> 
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
> 
> other attached packages:
> [1] vsn_3.12.0    Biobase_2.4.1 limma_2.18.3
> 
> loaded via a namespace (and not attached):
> [1] affy_1.22.1          affyio_1.12.0        grid_2.9.2           
> lattice_0.17-25      preprocessCore_1.6.0
>  >
> 
> ----------------------------------------------------------------------
> Dr. Tobias Straub ++4989218075439 Adolf-Butenandt-Institute, München D
> 
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: 
> http://news.gmane.org/gmane.science.biology.informatics.conductor

-- 
James W. MacDonald, M.S.
Biostatistician
Douglas Lab
University of Michigan
Department of Human Genetics
5912 Buhl
1241 E. Catherine St.
Ann Arbor MI 48109-5618
734-615-7826