Batch correction for RNA-seq didn't work with ComBat-seq tool

Question

[ChuYi]

Hello,

I was using ComBat-seq from the sva package to correct the batch effects for my RNA-seq samples, which were from two different batches but exactly the same platform.

Before doing this, I've checked my data for batch effects, here is the code I used and returned results:

library(DESeq2)
sample_info <- data.frame(sample = colnames(count_mat),
                            condition = factor(rep(c("pre_treat", "post_treat"),
                                                   times = 9)),
                            batch = factor(c(rep(1, 8), rep(2, 10))))
rownames(sample_info) <- sample_info$sample
dds <- DESeqDataSetFromMatrix(countData = count_mat,
                                colData = sample_info,
                                design= ~ condition+batch)
pca_dat <- plotPCA(DESeqTransform(dds), intgroup=c("condition", "batch"), 
                     returnData=TRUE)

I did the batch effects correction with ComBat-seq , and the input counts was the raw counts matrix, similarly, I also plotted the PCA:

gene_exp_adj <- ComBat_seq(counts = counts, batch = sample_info$batch, 
                             group = sample_info$condition)

but apparently, the results didn't change much:

I also performed same analysis on my lncRNA datasets, and the results of before/after batch correction were very similar. Does anyone know why? Or is there another tool I could try? Any advice would be appreciated!

[rpolicastro]

What do you mean when you say it didn't work. Was there an error, or was the batch effect not resolved?

[ChuYi]

The later case, the batch effects didn't seem to be resolved.

[4galaxy77]

Are you sure there’s any actual batch effects to begin with? I’m not an RNAseq expert but it doesn’t seem like it.

[ChuYi]

Maybe the results for my lncRNA datasets is more obvious, nothing seems to have changed: before batch correction

after batch correction

Answer 1

[ATpoint]

First of all you should follow the DESeq2 manual and use plotPCA correctly. You are giving it explicitely a DESeqTransform object (the manual does not suggest that -- it also makes no sense) and the axis limits of the PCA indicate that data are neither log-transformed - and based on the code probably not normalized as well. Please read the manual, do exactly what it does for the PCA and repeat the plots. My guess is that you are doing PCA on raw counts which makes no sense at all. So manual => DESeqDataSet() => vst() => plotPCA().

[ChuYi]

I've made a basic mistake, thank you so much⭐️⭐️

[ATpoint]

Is it fine now after applying these suggestions? If not feel free to ask for clarification.

[ChuYi]

Thank you for your enthusiastic answer! Well, after normalising the original data, I got the pca results as below:

gene datasets-before batch correction:

gene datasets-after batch correction:

Also, I tried removeBatchEffect function from limma package on gene datasets, but how to determine which one has better batch correction effect?

batch correction for gene datasets with limma:

But combat_seq didn't handle my lncRNA datasets well.

lncRNA datasets-before batch correction

lncRNA datasets-after batch correction

So I also tried removeBatchEffect function from limma package, it did seem to work:

batch correction for lncRNA datasets with limma:

Still not knowing the reason why lncRNA datasets have more obvious batch effect than mRNA datasets, since both of them were sequenced from the exact same source. And why combat_seq failed to correct batch effects for lncRNA datasets.