Hi bioconductorlist, I am trying to combine microarray data from different experiments(batches) where some of my batches fall completely in one condition (normal/cancer). For this I try to use ComBat in the sva package. After playing around with my data and ComBat it seemed to me that ComBat might introduce more signal in my data than there really is. The example given by ComBat-help (?ComBat) has a batch and condition setup much like my data, so I will use that to demonstrate. I will substitute the real data with random data ensuring no batch effect or condition effect. Then I will do a quick find- differentially-expressed-genes analysis with limma. Contradictory to what I would expect I end up with a list of genes. # start example code: # Starting off from ComBat-help (?ComBat) library(sva) library(bladderbatch) data(bladderdata) # Obtain phenotypic data pheno = pData(bladderEset) edata = exprs(bladderEset) batch = pheno$batch mod = model.matrix(~as.factor(cancer), data=pheno) # using dummy data with no effects instead of real data set.seed(100) edata[,] = rnorm(length(edata)) # Correct for batch using ComBat combat_edata = ComBat(dat=edata, batch=batch, mod=mod, par.prior=TRUE, prior.plots=FALSE) # overview of experimental layout with regards to conditions and batches # The "Biopsy" and "Normal" conditions seems to by most unbalanced. table(pheno[,3:4]) # Using the combat corrected dummy data to search for diff genes between "Biopsy"and "Normal" with limma library(limma) contrast = c("Biopsy-Normal") fac=factor(pheno$cancer) design = model.matrix(~0 + fac) colnames(design)=levels(fac) cont.matrix = makeContrasts ( contrasts=contrast, levels=design) # before ComBat fit <- lmFit(edata, design) fit2 = contrasts.fit(fit, cont.matrix) fit2 <- eBayes(fit2) tab = topTable(fit2, adjust='fdr', coef=contrast, number=999999) # top ten genes tab[1:10,] # no differetially exressed geens found print( paste("DUMMY DATA --- diff genes before batchnorm. Genes with FDR < 0.1: " , sum(tab$adj.P.Val<0.1), sep="")) # ComBat adjusted data fit <- lmFit(combat_edata, design) fit2 = contrasts.fit(fit, cont.matrix) fit2 <- eBayes(fit2) tab = topTable(fit2, adjust='fdr', coef=contrast, number=999999) # top ten genes tab[1:10,] # Many genes within 0.1 FDR. Looks similar to many microarray results. I would have reported this as a finding. print( paste("DUMMY DATA --- diff genes after batchnorm. Genes with FDR < 0.1: " , sum(tab$adj.P.Val<0.1), sep="")) # end example code I also tried (not shown) to use the real example data with permuted pheno data and run it through limma. The result was no diff-genes before ComBat and a lot after. And with fewer genes and an even more unbalanced setup I managed to make the standard hierarchical cluster heatmap before and after Combat with the samples clustering perfect according to condition afterwards(all with random start data). My interpretation of this is that starting with a unbalanced batch/condition setup like the one in the example (and in my data) and using ComBat will almost ensure that I end up with a list of significantly differentially expressed genes regardless of the real difference. This seems very wrong and I suspect I am missing something important. Please help me. Thanks for your help. Vegard Nygaard Bioinformatics Core Facility Department of Tumor Biology (Montebello) Institute for Cancer Research Oslo University Hospital Phone 0047 22781756 >sessionInfo() R version 3.0.2 (2013-09-25) Platform: x86_64-apple-darwin10.8.0 (64-bit) locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 attached base packages: [1] parallel stats graphics grDevices utils datasets methods base other attached packages: [1] limma_3.18.5 bladderbatch_1.0.7 Biobase_2.22.0 BiocGenerics_0.8.0 sva_3.8.0 mgcv_1.7-26 nlme_3.1-111 [8] corpcor_1.6.6 BiocInstaller_1.12.0 loaded via a namespace (and not attached): [1] grid_3.0.2 lattice_0.20-24 Matrix_1.0-14 tools_3.0.2