Why Batch effect removal with Combat-seq and DESeq2 give different results?

I’m trying to adjust batch effect using deseq2 limma::removeBatchEffect and also Combat-Seq. With limma version, I can clearly see the batch effect is removed, where I see control from Batch1 is together with the other 3 controls from Batch2.

###### Batch Correction with limma removeBatchEffect #######

dds <- DESeqDataSetFromMatrix(countData = data,
                              colData = coldata,
                              design = ~ Samplebatch + cond)

dds <- DESeq(dds)
dds

dds$batch <- as.numeric(dds$Samplebatch)
dds$cond
dds$batch

## vst after adding batch information
vsd_abc <- vst(dds, blind = FALSE)
vsd_abc$Samplebatch
vsd_abc$batch

design0 <- model.matrix(~cond, colData(vsd_abc))
assay(vsd_abc) <- limma::removeBatchEffect(assay(vsd_abc), vsd_abc$batch, design=design0) ### Batch Correction

And then I used vsd_abc for PCA visualization:

PCA after batch correction with limma::removeBatchEffect

But from this limma:removeBatchEffect function of DESeq2, I don’t get any batch corrected counts or TPM / FPKM expression data. I need this TPM for performing ssGSEA and various other analyses which are normalized for the length of the gene.

So, I decided to do the Batch correction with Combat-seq, through which you can get the Batch corrected counts and then I can convert that to TPM / FPKM and go with it. And the Combat-Seq was performed like below:


table(coldata$Samplebatch)

Batch1 Batch2 
    16      3 

table(coldata$cond)

Control  group1  group2  group3  group4  group5 
      4       3       3       3       3       3 

library(sva)
batch <- as.numeric(factor(coldata$Samplebatch))
cond2 <- as.numeric(factor(coldata$cond))
adjusted <- sva::ComBat_seq(as.matrix(data), batch=batch, group=cond2)

dds <- DESeqDataSetFromMatrix(countData = adjusted,
                              colData = coldata,
                              design = ~ Samplebatch + cond)

dds <- DESeq(dds)
dds$batch <- as.numeric(dds$Samplebatch)

## vst after adding batch information
vsd_abc <- vst(dds, blind = FALSE)

And then used vsd_abc for PCA that looked like below:

PCA after batch removal with Combat-seq

Here sample11 which is a Control from Batch1 is not clustered with other Controls from Batch2.

  1. Any idea what went wrong with Combat-seq batch removal?

  2. Is it possible to get TPM after batch effect removal using limma:removeBatchEffect?

Any help is appreciated. Thank you !!

Read more here: Source link