[Bioc-devel] Bug in rma: zero variance of probeset summaries across arrays

Thank you, Holger, Jim, Ben,

(1) this is indeed a peculiar feature of the medianpolish part of  
RMA. I apologize to Ben for using the harsh word "bug"  
inappropriately in this context. I was misled by observations that  
using different normalisation or background correction methods  
seemed to "fix" the problem, and had forgotten about those previous  
threads. Sorry.

I wonder whether anyone has done research on how to improve on  
medianpolish for small datasets and why people haven't adapted it.

(2) One further potential problem remains, and I hope that bringing  
this up might help to make sure that rma (which is obviously a core  
component of Bioconductor) is robust even under the most adverse  
circumstances (e.g. myself :). With the script below, I get a core  
dump. This is on a freshly checked out and built R (from today) and  
all packages then installed from scratch (sessionInfo as in previous  
mail below). It does not seem to happen on R-2.8.1. and Bioc release  
2.3., and even on R-dev, whether it happens seems to depend on the  
history of the R session. One hypothesis is that this could be  
related to the "destructive=TRUE" parameter of rma?

Can anyone reproduce this?

Best wishes
	Wolfgang


library("affy")
library("genefilter")
options(error=recover)

files = c("PP_T1_1.CEL", "PP_T1_2.CEL", "PP_T2_1.CEL", "PP_T2_2c.CEL")

#celdir = tempdir()
#for(f in files)
#  download.file(file.path("http://www.ebi.ac.uk/~huber/pub", f),
#                file.path(celdir, f))
celdir="."

a = ReadAffy(filenames=file.path(celdir, files))

## Define a CDF with one probeset per probe
testcdf = new.env(hash=TRUE)
for(i in seq_len(nrow(exprs(a))))
   assign(paste(i), matrix(i, ncol=1, nrow=1), envir=testcdf)
a at cdfName = "testcdf"

x = rma(a)


----------------------

Background correcting
Normalizing
Calculating Expression
Error in `colnames<-`(`*tmp*`, value = c("PP_T1_1.CEL",  
"PP_T1_2.CEL",  :
 attempt to set colnames on object with less than two dimensions

Enter a frame number, or 0 to exit

1: source("rma2.R")
2: eval.with.vis(ei, envir)
3: eval.with.vis(expr, envir, enclos)
4: rma(a)
5: `colnames<-`(`*tmp*`, value = c("PP_T1_1.CEL", "PP_T1_2.CEL",  
"PP_T2_1.CEL"

Selection:  4
Called from: eval(expr, envir, enclos)
Browse[1]>  ls()
[1] "*tmp*"       "background"  "bgversion"   "cols"         
"destructive"
[6] "exprs"       "ngenes"      "normalize"   "object"      "pNList"
[11] "rows"        "subset"      "verbose"
Browse[1]> str(exprs)

........ long dump of text ....


496: stop("attempt to set colnames on object with less than two  
dimensions")
497: `colnames<-`(`*tmp*`, value = c("PP_T1_1.CEL", "PP_T1_2.CEL",  
"PP_T2_1.CEL", "PP_T2_2c.CEL"))
498: rma(a)
499: eval.with.vis(expr, envir, enclos)
500: eval.with.vis(ei, envir)
501: source("rma2.R")

Possible actions:
1: abort (with core dump, if enabled)
2: normal R exit
3: exit R without saving workspace
4: exit R saving workspace
Selection:


James W. MacDonald wrote:

Hi Holger,
Holger Schwender wrote:

Hi Wolfgang,

the difference between the two outputs x1 and x2 are due to that in
x1 the probe intensities are background-corrected first and then
normalized, while in x2 they are first normalized and then
background-corrected.

The constant values are a feature of median polish: If median polish
is applied to four or six samples, then this can lead to constant
values for some of the probe sets. This has been mentioned years ago
by Jim (I couldn't find his original email on this topic in the
mailing list archive).

Actually this only occurs with a small, odd number of samples, not  
even. Well, it could happen for a small number of probesets if the  
n/2 and n/2 + 1 ordered values were identical, but it is much more  
common in the case of 3 or 5 samples.
Best,
Jim

Best, Holger


-------- Original-Nachricht --------

Datum: Tue, 03 Feb 2009 15:38:06 +0100 Von: Wolfgang Huber
<huber at ebi.ac.uk> An: Bioconductor Developers
<bioc-devel at stat.math.ethz.ch> Betreff: [Bioc-devel] Bug in rma:
zero variance of probeset summaries across    arrays

Hi,

I am getting puzzling results from the rma function in the affy
package, and furthermore they are different depending on whether
quantile normalisation is called within rma or explicitely. The
puzzling part is that some probesets, in a perfectly ordinary
experiment with 4 different arrays, are summarized with exactly
identical values across all arrays (i.e. zero variance); while the
underlying probe level data for these probesets is variable and
sits right in the middle of the dynamic range (i.e. not at any of
the extremes).

The code: ---------

library("affy") library("genefilter") library("preprocessCore")  
options(error=recover)

files = c("PP_T1_1.CEL", "PP_T1_2.CEL", "PP_T2_1.CEL",
"PP_T2_2c.CEL")

for(f in files) download.file(file.path("http://www.ebi.ac.uk/~huber/pub 
", f), file.path(tempdir(), f))

a = ReadAffy(filenames=file.path(tempdir(), files))

x1 = rma(a)

aq = a exprs(aq) = normalize.quantiles(exprs(aq)) x2 = rma(aq,
normalize=FALSE)

probeVars = function(x) { r = rowSds(exprs(x)) ps =
names(which(r==0)) cat(sprintf("%s: %d probes with zero
variance.\n", deparse(substitute(x)), length(ps))) return(ps) }

probeVars(a) ps = probeVars(x1) probeVars(x2)

print(exprs(a)[yeast2cdf[[ps[1]]][,1],])

-----------------------------------------------------------------------------


The output: ----------- Background correcting Normalizing  
Calculating Expression Background correcting Calculating Expression


a: 0 probes with zero variance. x1: 3 probes with zero variance.  
x2: 0 probes with zero variance.

PP_T1_1.CEL PP_T1_2.CEL PP_T2_1.CEL PP_T2_2c.CEL 35244          185
236         269          253 198153         226         314
320          213 48500          375         370         408
373 109094         125         157         201          175 144686
107         131         153          100 214110         181
296         307          234 47302          355         430
397          282 108656         117         130         150
118 147528         158         144         159          135 237691
208         244         253          169 229910          94
138         147          133

sessionInfo()

R version 2.9.0 Under development (unstable) (2009-02-03 r47829)  
x86_64-unknown-linux-gnu

locale:  
LC_CTYPE 
= 
it_IT 
.UTF 
-8 
;LC_NUMERIC 
= 
C 
;LC_TIME 
= 
it_IT 
.UTF 
-8 
;LC_COLLATE 
= 
it_IT 
.UTF 
-8 
;LC_MONETARY 
= 
C 
;LC_MESSAGES 
= 
it_IT 
.UTF 
-8 
;LC_PAPER 
= 
it_IT 
.UTF 
-8 
;LC_NAME 
= 
C 
;LC_ADDRESS 
=C;LC_TELEPHONE=C;LC_MEASUREMENT=it_IT.UTF-8;LC_IDENTIFICATION=C


attached base packages: [1] stats     graphics  grDevices datasets
utils     methods   base

other attached packages: [1] yeast2cdf_2.3.0
preprocessCore_1.5.3 genefilter_1.23.1 [4] affy_1.21.7
Biobase_2.3.10       fortunes_1.3-6

loaded via a namespace (and not attached): [1] affyio_1.11.3
annotate_1.21.3      AnnotationDbi_1.5.14 [4] DBI_0.2-4
RSQLite_0.7-1        splines_2.9.0 [7] survival_2.34-1
tools_2.9.0          xtable_1.5-4



Best wishes Wolfgang


PS While trying to dissect this, I also experimented with
constructing special CDF environments (e.g. one that had one
'probeset' for each probe), and I noted that repeatedly the C code
of rma caused core dumps. Perhaps a better response to unexpected
shapes of CDF environments would be to issue an error message. I
can provide details if interested.

---------------------------------------------------- Wolfgang
Huber, EMBL-EBI, http://www.ebi.ac.uk/huber

-- 

-- 
----------------------------------------------------
Wolfgang Huber  EMBL-EBI  http://www.ebi.ac.uk/huber