[BioC] possible problems with call.exprs for mas, mas-R, and justPlier

[Dick Beyer]

I'm comparing several BioC affy methods applied to the same data set.  I have noticed a discrepency between mas5(), call.exprs(method="mas"), and call.exprs(method="mas-R"), and a possible similar problem with justPlier.

Here is how I generated the data:

mydata.A.raw  <- ReadAffy()
mydata.B.raw  <- ReadAffy()

mydata.A.ce   <- call.exprs(mydata.A.raw,"mas5",sc=500)
mydata.B.ce   <- call.exprs(mydata.B.raw,"mas5",sc=500)
all.exprs.ce  <- rbind(exprs(mydata.A.ce),exprs(mydata.B.ce))

mydata.A.cer  <- call.exprs(mydata.A.raw,method="mas5-R",sc=500)
mydata.B.cer  <- call.exprs(mydata.B.raw,method="mas5-R",sc=500)
all.exprs.cer <- rbind(exprs(mydata.A.cer),exprs(mydata.B.cer))

mydata.A.mas  <- mas5(mydata.A.raw,sc=500)
mydata.B.mas  <- mas5(mydata.B.raw,sc=500)
all.exprs.mas <- rbind(log2(exprs(mydata.A.mas)),log2(exprs(mydata.B.mas)))

mydata.A.plr  <- justPlier(mydata.A.raw)
mydata.B.plr  <- justPlier(mydata.B.raw)
all.exprs.plr <- rbind(exprs(mydata.A.plr),exprs(mydata.B.plr))

all.exprs     <- NULL
all.exprs[1] <- list(all.exprs.ce)
all.exprs[2] <- list(all.exprs.cer)
all.exprs[3] <- list(all.exprs.mas)
all.exprs[4] <- list(all.exprs.plr)

design           <- model.matrix(~ -1+factor(c(rep(1,5),rep(2,5))))
colnames(design) <- c("group1","group2")

top.all  <- NULL
fit2.all <- NULL
for( i in 1:length(all.exprs) ){
  fit         <- lmFit(all.exprs[[i]], design)
  fit$Amean   <- rowMeans(all.exprs[[i]])
  cont.matrix <- makeContrasts(compare=group1-group2, levels=design)
  fit2        <- contrasts.fit(fit, cont.matrix)
  fit2        <- eBayes(fit2)
  fit2.all[i] <- list(fit2)
  top         <- topTable(fit2,coef=1,n=nrow(all.exprs[[i]]),
                   genelist=rownames(all.exprs[[i]]),
                   adjust.method="fdr",sort.by="P")
  o           <- order(top$Name)
  top.all[i]  <- list(top[o,])
}

Here are the summaries from these four calls:

> summary(fit2.all[[1]]$coef[,1])
     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-12.1900   0.1680   0.9118   1.0570   1.8030  10.9100 
> summary(fit2.all[[2]]$coef[,1])
     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-9.01800 -0.42810 -0.13710  0.02798  0.26990  7.95600 
> summary(fit2.all[[3]]$coef[,1])
     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-13.5100  -1.1550  -0.4108  -0.2666   0.4794   9.5860 
summary(fit2.all[[4]]$coef[,1])
     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-19.0700   0.2585   0.9192   1.1290   1.7710  19.2100

Both the call.exprs(method="mas") and justPlier have medians much different than zero.  The call to mas5() and call.exprs(method="mas-R") have medians closer to zero, but I would have thought they would be nearly the same.

I have tried all the various affy methods on this data set and the call.exprs(method="mas") and justPlier are different than the rest of the methods.

Perhaps I am doing something wrong in my use of call.exprs() or justPlier().

I am using HG-U133A/B arrays, and

base 2.0.1 
datasets 2.0.1 
utils 2.0.1 
grDevices 2.0.1 
graphics 2.0.1 
stats 2.0.1 
methods 2.0.1 
tools 2.0.1 
Biobase 1.5.0 
reposTools 1.5.1 
affy 1.5.8 
splines 2.0.1 
survival 2.16 
genefilter 1.5.0 
multtest 1.5.2 
siggenes 1.2.11 
qvalue 1.1 
limma 1.8.18 
matchprobes 1.0.12 
gcrma 1.1.1 
simpleaffy 2.08 
vsn 1.5.0 
hgu133acdf 1.4.3 
hgu133bcdf 1.4.3 
plier 0.03 
xtable 1.2-4


Thanks,
Dick
*******************************************************************************
Richard P. Beyer, Ph.D.	University of Washington
Tel.:(206) 616 7378	Env. & Occ. Health Sci. , Box 354695
Fax: (206) 685 4696	4225 Roosevelt Way NE, # 100
 			Seattle, WA 98105-6099
http://depts.washington.edu/ceeh/ServiceCores/FC5/FC5.html