[BioC] limma normexp background correction bug?
James W. MacDonald
jmacdon at med.umich.edu
Mon Oct 26 20:16:44 CET 2009
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
More information about the Bioconductor
mailing list