I recently started to use the backgroundCorrect function in limma
with the "normexp" method, as suggested in the user guide so as to
avoid the problem of zero or negative intensities. When analyzing
a set of about a dozen arrays, one of them ended up with NA M values
for every gene after the call to backgroundCorrect.
Looking at the code it transpires that backgroundCorrect calls
fit.normexp (when method="normexp"), and the code for fit.normexp
contains the following:
sigma <- sqrt(mean((f[f < mu] - mu)^2, na.rm = TRUE))
if (!is.infinite(sigma) || sigma < 1)
sigma <- 1
I.e. sigma is set to 1 unless it is +Inf, which I suspect might
not be what was intended. Removing the ! from is.infinite as a
test, and reanalysing the array which got set to all NAs, I do
get reasonable looking numbers back (not all NAs).
If this is a bug, I wonder what effect would it have had on other
data - would all data analyzed with backgroundCorrect in this way
be compromised?
Michael
Yes, this is a bug. The code was intended to be !is.finite(sigma)
rather than
!is.infinite(sigma). Thanks for finding this. Will be fixed in limma
version 1.8.14.
The bug affects only the starting values of fit.normexp(), not the
final parameter values. So
whenever the function converges, the bug will have had no final
effect. (I just tested the
function on a set of 34 arrays with and without the bug fix to confirm
this.)
What the bug will have done is to increase the frequency of
convergence failures. Convergence
failures lead to arrays with almost everything set to NA, so are
easily recognisable. limma has
been trapping and reporting convergence failures in fit.normexp()
since version 1.8.10. Hopefully
they will now occur less often.
Gordon
> Date: Wed, 15 Dec 2004 14:06:39 +1100
> From: michael_kirk@wmi.usyd.edu.au
> Subject: [BioC] Problem with limma backgroundCorrect function
> To: bioconductor@stat.math.ethz.ch
>
> I recently started to use the backgroundCorrect function in limma
> with the "normexp" method, as suggested in the user guide so as to
> avoid the problem of zero or negative intensities. When analyzing
> a set of about a dozen arrays, one of them ended up with NA M values
> for every gene after the call to backgroundCorrect.
>
> Looking at the code it transpires that backgroundCorrect calls
> fit.normexp (when method="normexp"), and the code for fit.normexp
> contains the following:
>
> sigma <- sqrt(mean((f[f < mu] - mu)^2, na.rm = TRUE))
> if (!is.infinite(sigma) || sigma < 1)
> sigma <- 1
>
> I.e. sigma is set to 1 unless it is +Inf, which I suspect might
> not be what was intended. Removing the ! from is.infinite as a
> test, and reanalysing the array which got set to all NAs, I do
> get reasonable looking numbers back (not all NAs).
>
> If this is a bug, I wonder what effect would it have had on other
> data - would all data analyzed with backgroundCorrect in this way
> be compromised?
>
> Michael
