[BioC] ComBat using SVA and bladderbatch package

[minyoung lee [guest]]

Dear users.

Hello? I have a question about ComBat results using SVA and bladderbatch package. The followings are the codes and the results I got.

library(sva)
library(pamr)
library(limma)

library(bladderbatch)
data(bladderdata)

pheno = pData(bladderEset)
edata = exprs(bladderEset)

mod = model.matrix(~as.factor(cancer), data=pheno)
mod0 = model.matrix(~1,data=pheno)

pValues = f.pvalue(edata,mod,mod0)
qValues = p.adjust(pValues,method="BH")
sum(qValues<=0.05)
[1] 15193
sum(qValues<=0.05)/length(qValues)
[1] 0.6818202

Note that nearly 70% of the genes are strongly differentially expressed at an FDR of less than
5% between groups. This number seems artificially high, even for a strong phenotype like cancer. This is the point that the authors emphasized in the bladderbatchTutorial. 

For ComBat,

batch = pheno$batch
mod = model.matrix(~as.factor(cancer), data=pheno)
combat_edata = ComBat(data=edata, batch=batch, mod=mod, numCovs=NULL, par.prior=TRUE, prior.plots=TRUE)

pValuesComBat = f.pvalue(combat_edata,mod,mod0)
qValuesComBat = p.adjust(pValuesComBat,method="BH")
sum(qValuesComBat<=0.05)
[1] 18300
sum(qValuesComBat<=0.05)/length(qValuesComBat)
[1] 0.8212539


After ComBat adjustment, 80% of the genes are differentially expressed at an FDR of less than 5% between groups. (The authors did not provide their results.) Is this result correct? 


 -- output of sessionInfo(): 

R version 2.14.2 (2012-02-29)
Platform: i386-pc-mingw32/i386 (32-bit)

locale:
[1] LC_COLLATE=Korean_Korea.949  LC_CTYPE=Korean_Korea.949   
[3] LC_MONETARY=Korean_Korea.949 LC_NUMERIC=C                
[5] LC_TIME=Korean_Korea.949    

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

other attached packages:
[1] bladderbatch_1.0.2 Biobase_2.14.0     limma_3.10.3      
[4] pamr_1.54          survival_2.37-4    cluster_1.14.3    
[7] sva_3.0.2          mgcv_1.7-22        corpcor_1.6.4     

loaded via a namespace (and not attached):
[1] grid_2.14.2     lattice_0.20-10 Matrix_1.0-5    nlme_3.1-108   


--
Sent via the guest posting facility at bioconductor.org.