limma::removeBatchEffect and using DESeq2 for RNA-seq differential expression

Question

[2405592M]

Hi everyone,

I have an RNA-seq experiment that I am analysing. I received cell culture replicates in batches. I used the limma::removeBatchEffect to remove any variance due to batch and re-did my PCA plots post batch correction, and it worked beautifully.

ddsObj$batch <- factor(rep(c("1", "2", "3", "4", "5")))
vsd <- varianceStabilizingTransformation(ddsObj)
plotPCA(vsd, "batch")
################################################
assay(vsd) <- limma::removeBatchEffect(assay(vsd), vsd$batch)    
plotPCA(vsd, intgroup = c("Condition"))

However, I used DESeq2 to do the differential expression and I've incorporated batch into my DESeq2 design:

design <- as.formula(~ batch + Condition)

My question is, even though I used Limma's remove batch effect to generate my lovely PCA plots (post DESeq2 analysis), would I be able to trust that DESeq2 removed the same variance generated by the batch effect that limma was so convincingly able to remove. Is my analysis downstream free from all the batch effects?

Have I done this the wrong way around? How can I confirm that the batch effect has been removed from my DESeq2 analysis or would I have to do my analysis in limma?

Answer 1

[leaodel]

You should 1)read the hundreds of posts about this and 2) dive into the docs to have a more thorough explanation, but yes: modeling a batch variable within the design is enough (and I must say more efficient than correcting) to have your batch effect unaccounted for in DE results. For visualization purposes, you'll need to remove the batch effect with something like limma.