Deleted:Why Batch effect removal with Combat-seq and DESeq2 give different results?
0
1
Entering edit mode
12 months ago
vk ▴ 10

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 !!

rnaseq combatseq batcheffect limma deseq2 • 1.1k views
ADD COMMENT
This thread is not open. No new answers may be added
Traffic: 1780 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6