Answers.Body: <p>For comparing the counts of different samples from DESeq2, Michael Love <a href="https://support.bioconductor.org/p/120764/#120843" rel="nofollow noreferrer">recommends</a> using the variance-stabilized transform. It'd be great if you could provide some specific code examples in your question, but without that here's something that should work with the DESeq2 workflow as mentioned in the package documentation:</p>
  
  <pre><code>## [assuming "dds &lt;- DESeq(dds)" has been run]
  dds.counts &lt;- assay(vst(dds, blind=FALSE));
  ## [This makes the scaled counts resemble more closely the actual read counts]
  dds.quantile99 &lt;- quantile(dds.counts[dds.counts &gt; min(dds.counts)], 0.99);
  dds.counts &lt;- ((dds.counts - min(dds.counts)) / (dds.quantile99-min(dds.counts))) * (dds.quantile99);
  </code></pre>
  
  <p>The matrix <em>dds.counts</em> should now represent a model-normalised approximation of the log of counts, with column names as the sample names, and row names as the gene names. The default behaviour of the VST changes zero counts to an arbitrary positive number, which is why I've modified the output so that count values of zero will be changed to 0 in log space (i.e. a count of 1). To plot results for a single gene, subset based on the row name:</p>
  
  <pre><code>gene.VST &lt;- sort(dds.counts["geneOfInterest",]);
  plot(gene.VST);
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-09-17T15:13:38.633Z
Answers.Text: Thanks! Sorry for not providing an example. I thought it was more about the practice of the community than a technical question so I omitted it. I use voom transformation followed by psva() in sva package but I assume a model-normalisation can also be done. The reason why I have this question is that using TPM shows a result closer to our expectation, while limma actually uses the voom transformation counts to detect differentially expressed genes. I am just not sure if I am misleading the reader if I use TPM just to show a better looking graph.
Answers.UserDisplayName: Kent
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-09-17T20:18:50.243Z
Answers.Text: Specific examples are helpful in bioinformatics because there are a huge number of solutions available that are context-dependent. It'd be great if you could add this information to your question to improve the quality of the answers you get.
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-09-17T21:22:01.717Z
Answers.Text: Hi! I have added the graphs that I get and the results. I am not sure if that's an example that would help answering the question, but that's basically the problem I am facing: will using TPM misleading/exaggerating the upregulation and trick people in interpreting the data wrong? And is it better to use voom since they are variance stabled and batch effect corrected so that reader can better interpret the data e.g. observe any trend in clinical data rather than randomness that might lead to the difference in gene A expression level within the group?
Answers.UserDisplayName: Kent
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-09-17T21:45:53.027Z
Answers.Text: It looks like you're doing a double normalization (i.e. using kallisto TPM values in addition to ``calcNormFactors(raw_data)``). I'm not familiar with voom, so don't know if voom is attempting additional correction despite the use of the ``normalize="none"`` option.
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-09-16T22:50:48.263Z
Answers.LastActivityDate: 2019-09-16T22:50:48.263Z
Answers.UserDisplayName: gringer
Answers.Body: <p>To add to what @gringer, when you do use TPM, the normalization done is for both library size and gene length. When you use rlog, the normalization is done via median normalization (<a href="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2864565/" rel="nofollow noreferrer">https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2864565/</a>). See below for a quick example about the rlog, and size factor normalization.</p>
  
  <pre><code>dds=makeExampleDESeqDataSet(betaSD=0.5,sizeFactors =runif(12))
  rld = rlog(dds)
  ## you can see the data is variance stablized, i.e lower counts are converted
  ### see the size factor calculated 
  rld$sizeFactor
  ###
  boxplot(assay(rld))
  ## look at unnormalized distributions
  boxplot(log2(counts(dds)+1))
  ## now we normalize and see
  dds = estimateSizeFactors(dds)
  # the size factors should be similar
  table(sizeFactors(dds) == rld$sizeFactor)
  ### here you can see the medians are similar but its a bit different from rlog
  boxplot(log2(counts(dds,normalize=TRUE)+1))
  ###
  </code></pre>
  
  <p>To answer your question specifically, the sizeFactor is unaware of any batch effect. It assumes median expression of all genes to be similar and normalize using that. If you think there is a batch effect, it has to be present across many genes, as opposed to one gene. I would check it with a PCA plot (<a href="https://www.rdocumentation.org/packages/DESeq2/versions/1.12.3/topics/plotPCA" rel="nofollow noreferrer">https://www.rdocumentation.org/packages/DESeq2/versions/1.12.3/topics/plotPCA</a>) and then test for differential expression by including the batch effect in the model.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-09-17T10:46:43.697Z
Answers.LastActivityDate: 2019-09-17T10:46:43.697Z
Answers.UserDisplayName: StupidWolf
Body: <p>I want to see the expression of a gene in a group of patient amongst the entire cohort using my RNA-Seq data. While I can do a differential expression analysis with limma or DESeq2, I want to see how much each sample from my cohort expresses the gene. </p>
  
  <p>The plan is to plot a waterfall plot (as defined in <a href="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4093310/" rel="nofollow noreferrer">this paper</a>). However, when I used TPM and rlog(CPM) they gave me very different looking graphs (Not just the shape, which is expected, but the ranking of expression). My concerns are, for TPM, the batch effect with two or three lots in one cohort, and, for rlog(CPM), the results don't match with prediction but is adjusted for batch effect. Which one would be a right or better choice?</p>
  
  <p>Upper one is with plot with. Lower one is plot with voom transformed counts (log(CPM)). I didn't include all cases since there are too many. This look more or less what it is like in the one including every case. Target is the condition group in which I am interested. Differential expression analysis definitely showed that it is upregulated.</p>
  
  <pre><code>                     logFC       CI.L       CI.R    AveExpr         t      P.Value    adj.P.Val         B
  GeneA  1.6576257  0.9143490  2.4009024  5.8713219  4.391796 1.646627e-05 0.0002581676  2.484413
  </code></pre>
  
  <p>Working steps
  TPM is basically from kallisto</p>
  
  <p>log(CPM) calculation</p>
  
  <pre><code>#Import data from tximport
  txi = tximport(files = read.hd5s,
                       type = "kallisto",
                       tx2gene = tx2gene,
                       countsFromAbundance = "lengthScaledTPM",
                       ignoreTxVersion = TRUE)
  raw_data = DGEList(txi$counts)
  # Filter low count genes
  keep &lt;- filterByExpr(raw_data)
  raw_data &lt;- raw_data[keep,,keep.lib.sizes=FALSE]
  # Normalise counts
  raw_data = calcNormFactors(raw_data)
  # log(CPM) by voom
  voom_data = voom(raw_data, design = NULL, normalize="none", plot = TRUE)
  # Batch effect correction
  voom_expr_psva = psva(as.matrix(voom_data$E),Batch)
  </code></pre>
  
  <p>Then I plot graph with TPM and log(CPM) respectively</p>
  
  <p><a href="https://i.sstatic.net/f1t1V.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/f1t1V.png" alt="Waterfall plot using TPM"></a>
  <a href="https://i.sstatic.net/jxRIp.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/jxRIp.png" alt="Waterfall plot using voom"></a></p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-09-17T21:40:23.883Z
Comments.Text: in addition to TPM, I think kallisto also produces count-like data as output, which would be better for downstream analysis (presumably the default when tximport is used in 'kallisto' mode). The count-level data gives a better idea of uncertainty based on shot noise (i.e. the likelihood of missing a gene because its expression levels are low).
Comments.UserDisplayName: gringer
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-09-17T22:01:27.097Z
Comments.Text: I think there might be some misunderstanding. I actually used the count data as shown in this line: raw_data = DGEList(txi$counts), which is to generate a DGEList object from the count-like data you mentioned.
Comments.UserDisplayName: Kent
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-09-17T23:04:59.150Z
Comments.Text: You're using the `countsFromAbundance="lengthScaledTPM"` argument on tximport, which looks suspiciously like a form of normalisation to me ("Whether to generate estimated counts using abundance estimates"). The import default for that, as used in the [kallisto import example](https://bioconductor.org/packages/release/bioc/vignettes/tximport/inst/doc/tximport.html#kallisto) is `countsFromAbundance="no"`.
Comments.UserDisplayName: gringer
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-09-18T13:16:21.210Z
Comments.Text: Thanks for pointing that out. I didn't notice since it was something shown to me by a bioinformatician when I first started learning R. Yes it looks like the reads are scaled to the library size and average transcript length as show in this [tutorial](https://bioconductor.org/packages/devel/bioc/vignettes/tximport/inst/doc/tximport.html). I guess `raw_data = calcNormFactors(raw_data)` was intended to renormalise to library size after filtering due to `keep.lib.sizes = FALSE`. Still trying to work this around my head.
Comments.UserDisplayName: Kent
Tags: <rna-seq><limma><batch-effects><tpm>
Title: TPM or rlog(CPM) for comparing expression?
ThreadText: # TPM or rlog(CPM) for comparing expression?
  
  (Asked by: username_0 on 2019-09-16)
  
  I want to see the expression of a gene in a group of patient amongst the entire cohort using my RNA-Seq data. While I can do a differential expression analysis with limma or DESeq2, I want to see how much each sample from my cohort expresses the gene.
  
  The plan is to plot a waterfall plot (as defined in [this paper](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4093310/)). However, when I used TPM and rlog(CPM) they gave me very different looking graphs (Not just the shape, which is expected, but the ranking of expression). My concerns are, for TPM, the batch effect with two or three lots in one cohort, and, for rlog(CPM), the results don't match with prediction but is adjusted for batch effect. Which one would be a right or better choice?
  
  Upper one is with plot with. Lower one is plot with voom transformed counts (log(CPM)). I didn't include all cases since there are too many. This look more or less what it is like in the one including every case. Target is the condition group in which I am interested. Differential expression analysis definitely showed that it is upregulated.
  
  ```
                       logFC       CI.L       CI.R    AveExpr         t      P.Value    adj.P.Val         B
  GeneA  1.6576257  0.9143490  2.4009024  5.8713219  4.391796 1.646627e-05 0.0002581676  2.484413
  
  ```
  
  Working steps
  TPM is basically from kallisto
  
  log(CPM) calculation
  
  ```
  #Import data from tximport
  txi = tximport(files = read.hd5s,
                       type = "kallisto",
                       tx2gene = tx2gene,
                       countsFromAbundance = "lengthScaledTPM",
                       ignoreTxVersion = TRUE)
  raw_data = DGEList(txi$counts)
  # Filter low count genes
  keep <- filterByExpr(raw_data)
  raw_data <- raw_data[keep,,keep.lib.sizes=FALSE]
  # Normalise counts
  raw_data = calcNormFactors(raw_data)
  # log(CPM) by voom
  voom_data = voom(raw_data, design = NULL, normalize="none", plot = TRUE)
  # Batch effect correction
  voom_expr_psva = psva(as.matrix(voom_data$E),Batch)
  
  ```
  
  Then I plot graph with TPM and log(CPM) respectively
  
  [![Waterfall plot using TPM](https://i.sstatic.net/f1t1V.png)](https://i.sstatic.net/f1t1V.png)
  [![Waterfall plot using voom](https://i.sstatic.net/jxRIp.png)](https://i.sstatic.net/jxRIp.png)
  
  **Tags:** rna-seq, limma, batch-effects, tpm
  
  > **Comments:**
  > **username_1**: in addition to TPM, I think kallisto also produces count-like data as output, which would be better for downstream analysis (presumably the default when tximport is used in 'kallisto' mode). The count-level data gives a better idea of uncertainty based on shot noise (i.e. the likelihood of missing a gene because its expression levels are low).  
  > **username_0**: I think there might be some misunderstanding. I actually used the count data as shown in this line: raw_data = DGEList(txi$counts), which is to generate a DGEList object from the count-like data you mentioned.  
  > **username_1**: You're using the `countsFromAbundance="lengthScaledTPM"` argument on tximport, which looks suspiciously like a form of normalisation to me ("Whether to generate estimated counts using abundance estimates"). The import default for that, as used in the [kallisto import example](https://bioconductor.org/packages/release/bioc/vignettes/tximport/inst/doc/tximport.html#kallisto) is `countsFromAbundance="no"`.  
  > **username_0**: Thanks for pointing that out. I didn't notice since it was something shown to me by a bioinformatician when I first started learning R. Yes it looks like the reads are scaled to the library size and average transcript length as show in this [tutorial](https://bioconductor.org/packages/devel/bioc/vignettes/tximport/inst/doc/tximport.html). I guess `raw_data = calcNormFactors(raw_data)` was intended to renormalise to library size after filtering due to `keep.lib.sizes = FALSE`. Still trying to work this around my head.  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_1 on 2019-09-16. Score: 3
  
  For comparing the counts of different samples from DESeq2, Michael Love [recommends](https://support.bioconductor.org/p/120764/#120843) using the variance-stabilized transform. It'd be great if you could provide some specific code examples in your question, but without that here's something that should work with the DESeq2 workflow as mentioned in the package documentation:
  
  ```
  ## [assuming "dds <- DESeq(dds)" has been run]
  dds.counts <- assay(vst(dds, blind=FALSE));
  ## [This makes the scaled counts resemble more closely the actual read counts]
  dds.quantile99 <- quantile(dds.counts[dds.counts > min(dds.counts)], 0.99);
  dds.counts <- ((dds.counts - min(dds.counts)) / (dds.quantile99-min(dds.counts))) * (dds.quantile99);
  
  ```
  
  The matrix *dds.counts* should now represent a model-normalised approximation of the log of counts, with column names as the sample names, and row names as the gene names. The default behaviour of the VST changes zero counts to an arbitrary positive number, which is why I've modified the output so that count values of zero will be changed to 0 in log space (i.e. a count of 1). To plot results for a single gene, subset based on the row name:
  
  ```
  gene.VST <- sort(dds.counts["geneOfInterest",]);
  plot(gene.VST);
  
  ```
  
  > **Comments:**
  > **username_0**: Thanks! Sorry for not providing an example. I thought it was more about the practice of the community than a technical question so I omitted it. I use voom transformation followed by psva() in sva package but I assume a model-normalisation can also be done. The reason why I have this question is that using TPM shows a result closer to our expectation, while limma actually uses the voom transformation counts to detect differentially expressed genes. I am just not sure if I am misleading the reader if I use TPM just to show a better looking graph.  
  > **username_1**: Specific examples are helpful in bioinformatics because there are a huge number of solutions available that are context-dependent. It'd be great if you could add this information to your question to improve the quality of the answers you get.  
  > **username_0**: Hi! I have added the graphs that I get and the results. I am not sure if that's an example that would help answering the question, but that's basically the problem I am facing: will using TPM misleading/exaggerating the upregulation and trick people in interpreting the data wrong? And is it better to use voom since they are variance stabled and batch effect corrected so that reader can better interpret the data e.g. observe any trend in clinical data rather than randomness that might lead to the difference in gene A expression level within the group?  
  > **username_1**: It looks like you're doing a double normalization (i.e. using kallisto TPM values in addition to ``calcNormFactors(raw_data)``). I'm not familiar with voom, so don't know if voom is attempting additional correction despite the use of the ``normalize="none"`` option.  
  
  ---
  
  ### Answer by username_2 on 2019-09-17. Score: 2
  
  To add to what @username_1, when you do use TPM, the normalization done is for both library size and gene length. When you use rlog, the normalization is done via median normalization (<https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2864565/>). See below for a quick example about the rlog, and size factor normalization.
  
  ```
  dds=makeExampleDESeqDataSet(betaSD=0.5,sizeFactors =runif(12))
  rld = rlog(dds)
  ## you can see the data is variance stablized, i.e lower counts are converted
  ### see the size factor calculated 
  rld$sizeFactor
  ###
  boxplot(assay(rld))
  ## look at unnormalized distributions
  boxplot(log2(counts(dds)+1))
  ## now we normalize and see
  dds = estimateSizeFactors(dds)
  # the size factors should be similar
  table(sizeFactors(dds) == rld$sizeFactor)
  ### here you can see the medians are similar but its a bit different from rlog
  boxplot(log2(counts(dds,normalize=TRUE)+1))
  ###
  
  ```
  
  To answer your question specifically, the sizeFactor is unaware of any batch effect. It assumes median expression of all genes to be similar and normalize using that. If you think there is a batch effect, it has to be present across many genes, as opposed to one gene. I would check it with a PCA plot (<https://www.rdocumentation.org/packages/DESeq2/versions/1.12.3/topics/plotPCA>) and then test for differential expression by including the batch effect in the model.
  
  
  ---
  
  