Body: <p>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 <code>PCA</code> plot and seeing the batch effect in PC1, I used <code>removeBatchEffect</code> 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 <code>model.matrix</code>. So, I did that, ran <code>limma/voom</code> 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 <code>removeBatchEffect</code> on the count data from DEGs and then use the transformed dataset for heatmaps or there is another way to fix this?</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-08-15T13:01:44.627Z
Comments.Text: Maybe you could try more state-of-the-art of methods like "comBat" from "sva" package to remove the batch effects.
Comments.UserDisplayName: Phoenix Mu
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-08-15T13:14:22.187Z
Comments.Text: 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
Comments.UserDisplayName: Reza Rezaei
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-08-15T13:16:34.350Z
Comments.Text: Read [this](https://www.biostars.org/p/266507/) for more information.
Comments.UserDisplayName: Reza Rezaei
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-08-16T13:56:22.997Z
Comments.Text: 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.
Comments.UserDisplayName: Phoenix Mu
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-08-16T15:53:12.660Z
Comments.Text: 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?
Comments.UserDisplayName: Reza Rezaei
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-08-16T17:32:59.080Z
Comments.Text: If your data is log2 scaled, are there negative values before removing any batch effects?
Comments.UserDisplayName: Phoenix Mu
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-08-16T17:40:15.820Z
Comments.Text: I just checked and `Yes` there are a few negative values before removing batch effects.
Comments.UserDisplayName: Reza Rezaei
Tags: <rna-seq><edger><limma><heatmap><batch-effects>
Title: Removing Batch Effect in Heatmaps after Differential Gene Expression Analysis
ThreadText: # Removing Batch Effect in Heatmaps after Differential Gene Expression Analysis
  
  (Asked by: username_0 on 2020-08-15)
  
  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?
  
  **Tags:** rna-seq, edger, limma, heatmap, batch-effects
  
  > **Comments:**
  > **username_1**: Maybe you could try more state-of-the-art of methods like "comBat" from "sva" package to remove the batch effects.  
  > **username_0**: 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  
  > **username_0**: Read [this](https://www.biostars.org/p/266507/) for more information.  
  > **username_1**: 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.  
  > **username_0**: 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?  
  > **username_1**: If your data is log2 scaled, are there negative values before removing any batch effects?  
  > **username_0**: I just checked and `Yes` there are a few negative values before removing batch effects.  
  
  ---
  
  ## 0 answers
  
  