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:
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:
Here sample11
which is a Control from Batch1 is not clustered with other Controls from Batch2.
-
Any idea what went wrong with
Combat-seq
batch removal? -
Is it possible to get TPM after batch effect removal using
limma:removeBatchEffect
?
Any help is appreciated. Thank you !!
Read more here: Source link