Removing Batch Effect in Heatmaps after Differential Gene Expression Analysis
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.