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?
1 Answer
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-seqfor 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-seqand 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-seqon my raw data and then did the normalization steps and compared its PCA plot with the corresponding plot forcomBat. It not only didn't improve the distance between replicates but, not significantly,comBatseems 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$– user3051Aug 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$– AhdeeAug 13, 2021 at 18:31
removeBatchEffectfromlimmaovercomBatfromsva. The gist of it was that we should not use tools likecomBatbecause they are biased towards removing all batch effects. Moreover,comBatand alsoremoveBatchEffectgenerate 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 useremoveBatchEffectfor visualizations but just use batch inmodel.matrixfor differential gene expression analysis $\endgroup$comBatand alsoremoveBatchEffect. 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$