Answers.Body: <p>There are a number of methods. If you're doing DE, you have ComBatSeq, SVAseq, RUVseq, BUSseq. You could also try Z-score normalization (<span class="math-container">$(x - \overline{x}) / \sigma$</span>) or even quantile transformation. For the latter two, make sure you work on each batch individually, not the whole dataset at once. See more on that <a href="https://www.nature.com/articles/s41598-020-72664-6" rel="nofollow noreferrer">here</a>.</p>
  <p>To visually compare if transformation works, plotting PCA/UMAP/t-SNE on raw and transformed data can perhaps be of some insight.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-04-21T08:44:08.253Z
Answers.Text: Thanks Roman for the response. Very useful. Do you recommend to trim the reads to 125bp for the NovaSeq to compare with the other batches from HiSeq? Or do you suggest to include read length as a factor during DE analysis? What do you suggest.
Answers.UserDisplayName: Praveen
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-04-22T07:13:21.133Z
Answers.Text: This sounds like another question. :) Feel free to post it as such.
Answers.UserDisplayName: Roman Luštrik
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-04-22T08:29:32.893Z
Answers.Text: Okay, let me do that. Don't know if this kind of experiment has been done previously (happens data is used from multiple sources for integrated analysis)
Answers.UserDisplayName: Praveen
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-04-21T05:23:08.627Z
Answers.LastActivityDate: 2021-04-21T05:23:08.627Z
Answers.UserDisplayName: Roman Luštrik
Answers.Body: <p>I agree with <a href="https://bioinformatics.stackexchange.com/users/2388/swbarnes2">swbarnes</a> that you don't need to explicitly model for different read length, but you need to normalize for different production of sequencing runs.
  Among the several available packages, <a href="https://bioconductor.org/packages/release/bioc/html/DESeq2.html" rel="nofollow noreferrer">DEseq2</a> is one of the most widely used.
  The normalization method computes the geometric mean of expression of each gene across samples; then, expression of the gene in each sample is divided by the geometric mean. Once this has been don for all genes, in each sample it computes the median of values, and this is the &quot;size factor&quot;. Ok, it looks convoluted, but it is indeed easy, and it is implemented in DESeq2 and in other software packages.
  You can find methodological details <a href="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4302049/" rel="nofollow noreferrer">here</a>.
  Finally, DESeq2 also offer the possibility to account for additional experimental details (e.g. read length) and this should reduce the risk of spurious results.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-04-26T10:47:59.787Z
Answers.Text: Thanks for your response, Fabio. If I use ether median-to-ratios (geometric mean in DEseq2) or upper-quartile normalization, I think it should work good despite different sequencing machines? Also do you have any articles where the read length is used as a factor in the design for DE analysis? How robust is this method compared to trimming of reads as trimming on large samples is a huge additional task.
Answers.UserDisplayName: Praveen
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-04-24T13:05:52.923Z
Answers.LastActivityDate: 2021-04-24T13:12:05.420Z
Answers.LastEditDate: 2021-04-24T13:12:05.420Z
Answers.UserDisplayName: Fabio Marroni
Answers.Body: <p>Normalization and batch-effect correction are two different things.</p>
  <p>You always have to do normalization of RNAseq data in order to be able to compare the runs and do DE analysis, for instance. Plenty of methods and you have to search for them as normalization methods indeed.</p>
  <p>What you are asking and concerned about is &quot;batch-effects&quot; in your data as a function of having used two different Illumina machines. This is a valid concern, but before you try to correct for it or mitigate it, you first have to assess whether it is actually there. Maybe is not at all. To figure out if you have batch-effects you have to first do data exploration, this is why the use of PCA/UMAP/t-SNE has been already suggested. In this plots you will be trying to see whether your data cluster together by machine type instead of by the biology it represents. Read about how to assess the presence of batch effects.</p>
  <p>If you do observe it or anyway want to correct for its presence, then indeed you can try one of the methods also already mentioned: ComBatSeq, SVAseq, RUVseq, BUSseq.</p>
  <p>Rather than the two different machines, the batch-effect might actually be due to different people doing the lib preparation, experimental settings if experiments were done at different moments, etc.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-04-21T09:18:10.810Z
Answers.LastActivityDate: 2021-04-23T14:10:18.647Z
Answers.LastEditDate: 2021-04-23T14:10:18.647Z
Answers.UserDisplayName: JRodrigoF
Answers.Body: <p>I don't think RNASEq is that sensitive to running different things on different instruments.  Differing read lengths will change things.  You could always trim the reads down so they all match.  Or include read length as a factor in the design when looking for DE genes.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-04-20T17:37:11.240Z
Answers.LastActivityDate: 2021-04-20T17:37:11.240Z
Answers.UserDisplayName: swbarnes2
Body: <p>I have RNA samples sequenced with TruSeq Stranded Total RNA kit protocol in Illumina HiSeq (2x125bp) and NovaSeq platforms (2x150bp) - almost 100 samples altogether. I have to use the samples data for expression analysis. Question is how do we normalize samples from different runs and of read length? What is best way (with available packages) we can normalize the count (to overcome the bias) and do a DE analysis? ANy suggestion would be helpful.</p>
  
Tags: <rna-seq><differential-expression><gene-expression><normalization><featurecounts>
Title: Normalize RNA seq data from multiple runs for expression analysis
ThreadText: # Normalize RNA seq data from multiple runs for expression analysis
  
  (Asked by: username_0 on 2021-04-20)
  
  I have RNA samples sequenced with TruSeq Stranded Total RNA kit protocol in Illumina HiSeq (2x125bp) and NovaSeq platforms (2x150bp) - almost 100 samples altogether. I have to use the samples data for expression analysis. Question is how do we normalize samples from different runs and of read length? What is best way (with available packages) we can normalize the count (to overcome the bias) and do a DE analysis? ANy suggestion would be helpful.
  
  **Tags:** rna-seq, differential-expression, gene-expression, normalization, featurecounts
  
  
  ---
  
  ## 4 answers
  
  ### Answer by username_1 on 2021-04-21. Score: 1
  
  There are a number of methods. If you're doing DE, you have ComBatSeq, SVAseq, RUVseq, BUSseq. You could also try Z-score normalization ($(x - \overline{x}) / \sigma$) or even quantile transformation. For the latter two, make sure you work on each batch individually, not the whole dataset at once. See more on that [here](https://www.nature.com/articles/s41598-020-72664-6).
  
  To visually compare if transformation works, plotting PCA/UMAP/t-SNE on raw and transformed data can perhaps be of some insight.
  
  > **Comments:**
  > **username_0**: Thanks username_1 for the response. Very useful. Do you recommend to trim the reads to 125bp for the NovaSeq to compare with the other batches from HiSeq? Or do you suggest to include read length as a factor during DE analysis? What do you suggest.  
  > **username_1**: This sounds like another question. :) Feel free to post it as such.  
  > **username_0**: Okay, let me do that. Don't know if this kind of experiment has been done previously (happens data is used from multiple sources for integrated analysis)  
  
  ---
  
  ### Answer by username_2 on 2021-04-24. Score: 0
  
  I agree with [swbarnes](https://bioinformatics.stackexchange.com/users/2388/username_4) that you don't need to explicitly model for different read length, but you need to normalize for different production of sequencing runs.
  Among the several available packages, [DEseq2](https://bioconductor.org/packages/release/bioc/html/DESeq2.html) is one of the most widely used.
  The normalization method computes the geometric mean of expression of each gene across samples; then, expression of the gene in each sample is divided by the geometric mean. Once this has been don for all genes, in each sample it computes the median of values, and this is the "size factor". Ok, it looks convoluted, but it is indeed easy, and it is implemented in DESeq2 and in other software packages.
  You can find methodological details [here](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4302049/).
  Finally, DESeq2 also offer the possibility to account for additional experimental details (e.g. read length) and this should reduce the risk of spurious results.
  
  > **Comments:**
  > **username_0**: Thanks for your response, username_2. If I use ether median-to-ratios (geometric mean in DEseq2) or upper-quartile normalization, I think it should work good despite different sequencing machines? Also do you have any articles where the read length is used as a factor in the design for DE analysis? How robust is this method compared to trimming of reads as trimming on large samples is a huge additional task.  
  
  ---
  
  ### Answer by username_3 on 2021-04-21. Score: 0
  
  Normalization and batch-effect correction are two different things.
  
  You always have to do normalization of RNAseq data in order to be able to compare the runs and do DE analysis, for instance. Plenty of methods and you have to search for them as normalization methods indeed.
  
  What you are asking and concerned about is "batch-effects" in your data as a function of having used two different Illumina machines. This is a valid concern, but before you try to correct for it or mitigate it, you first have to assess whether it is actually there. Maybe is not at all. To figure out if you have batch-effects you have to first do data exploration, this is why the use of PCA/UMAP/t-SNE has been already suggested. In this plots you will be trying to see whether your data cluster together by machine type instead of by the biology it represents. Read about how to assess the presence of batch effects.
  
  If you do observe it or anyway want to correct for its presence, then indeed you can try one of the methods also already mentioned: ComBatSeq, SVAseq, RUVseq, BUSseq.
  
  Rather than the two different machines, the batch-effect might actually be due to different people doing the lib preparation, experimental settings if experiments were done at different moments, etc.
  
  
  ---
  
  ### Answer by username_4 on 2021-04-20. Score: 0
  
  I don't think RNASEq is that sensitive to running different things on different instruments. Differing read lengths will change things. You could always trim the reads down so they all match. Or include read length as a factor in the design when looking for DE genes.
  
  
  ---
  
  