I'm working on a dataset in which the first replicate of each group is one batch and the second replicate is in a second batch. After checking the PCA plot and seeing the batch effect in PC1, I used removeBatchEffect function from limma to subtract the batch effect from my count data. Then, using PCA gives me a plot that doesn't seem to have any apparent batch effect left! However, it is recommended not to use batch effect correction for differential gene analysis but use the batch variable along with the group variable in constructing the model.matrix. So, I did that, ran limma/voom on the normalized counts, and extracted differentially expressed genes. However, when I'm trying to make a heatmap from the DEGs, I still see that the samples from different batches are clustered separately, instead of seeing the clustering of replicates of the same sample. So, my question is whether I should use removeBatchEffect on the count data from DEGs and then use the transformed dataset for heatmaps or there is another way to fix this?

  • $\begingroup$ Maybe you could try more state-of-the-art of methods like "comBat" from "sva" package to remove the batch effects. $\endgroup$
    – Phoenix Mu
    Aug 15, 2020 at 13:01
  • $\begingroup$ I read a few posts in biostar that confirmed my choice of removeBatchEffect from limma over comBat from sva. The gist of it was that we should not use tools like comBat because they are biased towards removing all batch effects. Moreover, comBat and also removeBatchEffect generate negative sign results which can't be processed by voom or edgeR and trying to remove these negative results is another level of manipulation in the dataset. So, it is better to use removeBatchEffect for visualizations but just use batch in model.matrix for differential gene expression analysis $\endgroup$ Aug 15, 2020 at 13:14
  • $\begingroup$ Read this for more information. $\endgroup$ Aug 15, 2020 at 13:16
  • 1
    $\begingroup$ It is true that comBat introduces negative values. This is because comBat should be used on a normalized matrix. If all values in the matrix is positive, I am assuming it is a count matrix (which a newer version of sva takes care of) or a TPM matrix (which should be normalized first.) I think maybe you could run comBat after voom if you decide to use that. But of course using limma in a right way could work as well. $\endgroup$
    – Phoenix Mu
    Aug 16, 2020 at 13:56
  • $\begingroup$ Could you please clarify your comment a bit. my data is already log2 scaled, filtered for lowly expressed genes and also TMM normalized, but I still get negative values with comBat and also removeBatchEffect. However, I'm using batch as a covariat in my model.matrix and do the DEG analysis. You meant I can use comBat on DEGs for visualizations? $\endgroup$ Aug 16, 2020 at 15:53

1 Answer 1


It is indeed true that for the DE analysis one should include batch into the formula to avoid changing the original counts. Still, for everything else such as plotting heatmaps use of removeBatchEffects is perfectly fine and (at least to me) a standard and well-accepted procedure. It essentially does not matter what you use to correct for the batch effect for the counts you use downstream. Results will probably be similar. Combat-Seq from the sva package is a recent adaptation of ComBat specifically for RNA-seq which (from what I understand) better deals with the integer-count nature of the data. This operates on raw counts and avoids the infamous negative values that happen at times with both limma and Combat. After applying ComBat-Seq on your raw counts you can normalize them as usual with edgeR (or any tool you like) and then make the heatmaps. See https://github.com/zhangyuqing/ComBat-seq. The DE results still should come from the normal DE pipeline with batch as covariate as discussed above.

  • $\begingroup$ Thank you for your answer. Just one confusing thing about it. I first normalize, filter, and log2 scale and then do the EdgeR analysis with batch as covariat. Then, should I use the normalized.filtered.log2 as the output for ComBat-seq for module identification using heatmaps or should I just select differentially expressed genes from my raw data and use that as input for heatmaps and clustering? $\endgroup$ Aug 16, 2020 at 17:03
  • $\begingroup$ Or maybe it is better to first use ComBat-seq and eliminate the batch effect and then go through my normal workflow? Is this your recommended approach? $\endgroup$ Aug 16, 2020 at 17:06
  • $\begingroup$ I just checked combat-seq on my raw data and then did the normalization steps and compared its PCA plot with the corresponding plot for comBat. It not only didn't improve the distance between replicates but, not significantly, comBat seems to perform better!! $\endgroup$ Aug 16, 2020 at 17:47
  • 1
    $\begingroup$ As I said above, do your standard edgeR analysis with batch as covariate to get DE genes. Then, as a second step to get batch-corrected counts do ComBat-Seq on the raw count matix, normalize and log2 it and use this for heatmaps (and every other downstream analysis). Without full code and plots I cannot comment on what you describe. In my hands ComBat-seq worked very well and removes batch effects effectively. $\endgroup$
    – user3051
    Aug 16, 2020 at 18:07
  • $\begingroup$ @ATpoint thanks I agree with your points. Do you think its fair to run the standard DE with the batch as covariate, but when generating the post analysis plot, to merge the corrected dataframe ( say in log cpm) back with the statistical output for plotting? $\endgroup$
    – Ahdee
    Aug 13, 2021 at 18:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.