Hello everyone,

I am running differential expression analysis between patients and healthy control

In the expression dataset the biological variable of interest and batch are partially confounded. Even though batch effect is not particularly strong on Principal components 1 versus 2, I would prefer to run some batch adjustment. However, when I run ComBat WITHOUT including my disease indicator in the model (according to the updated sva vignette)  the number of differentially expressed genes is 5 vs 28 without running combat.

I am concerned that batch correction have removed actual biological differences. What approach should I use to feel
more comfortable about a possible batch effect inflating my results?

Thank you

Adam Jaros

To illustrate I added 1) batches in my actual dataset 2) what I want to do - illustrated on bladderbatch data

1)batches in my actuall dataset

 2) To show the difference

library(sva)
library(oligo)
library(bladderbatch)
library(limma)

# load data
data(bladderdata)
p<-pData(bladderEset)


# create subset similar to what I have
use<-sampleNames(bladderEset)[p$outcome=='Normal' | p$outcome=='sTCC-CIS']
subEset<-bladderEset[,use]
p<-pData(subEset)
xtabs(~batch+cancer,p)

# create model for combat and run it
mod<-model.matrix(~1,data=p)
batch<-p$batch
combat<-subEset # duplicate subEset
exprs(combat)<-ComBat(dat=exprs(subEset),batch=batch,mod=mod)

groups<-p$outcome=='sTCC-CIS' #define groups for comparison
par(mfrow=c(1,2)) # graphics

# create a volcano plot for each of them
for (eset in list(combat,subEset)){
        p<-pData(eset)
        e<-exprs(eset)
        d<-rowMeans(e[,groups],na.rm=T)-rowMeans(e[,!groups],na.rm=T)
        
        # do row ttests
        tt<-rowttests(e,factor(groups))
        lod<--log10(tt[['p.value']])
        
        ## volcano plots
        smoothScatter(d,lod,xlim=c(-3,3),ylim=c(0,10))  
        abline(h=2,v=c(-1,1))
        
        
      
}

3) bladderbatch subset

     cancer
batch Biopsy Cancer Normal
    2      0     13      4
    3      0      0      4
    5      0      3      0

Hi Adam,

Sorry for the delay in response, I've been traveling. One thing to keep in mind is that when batch and the variable you care about are correlated, adjusting for batch effects will very frequently reduce the number of genes called differentially expressed. The danger of skipping batch adjustment is that some of the genes called differentially expressed may only be differentially expressed due to batch effects. 

It is not surprising to see the # of significant genes go from 28 to 5 when adjusting for batch effects. One thing you can do is plot the expression levels for each of the 28 genes colored by batch and colored by the variable you found DE and try to determine if it is possible to distinguish batch effects from effects you care about - but this is usually really dicey and rarely more informative than the Combat analysis in my experience. 

Jeff

Hi Adam,

Sorry for the delay in response, I've been traveling. One thing to keep in mind is that when batch and the variable you care about are correlated, adjusting for batch effects will very frequently reduce the number of genes called differentially expressed. The danger of skipping batch adjustment is that some of the genes called differentially expressed may only be differentially expressed due to batch effects. 

It is not surprising to see the # of significant genes go from 28 to 5 when adjusting for batch effects. One thing you can do is plot the expression levels for each of the 28 genes colored by batch and colored by the variable you found DE and try to determine if it is possible to distinguish batch effects from effects you care about - but this is usually really dicey and rarely more informative than the Combat analysis in my experience. 

Jeff