<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>

CC BY-SA 4.0
2019-09-17T15:13:38.633Z
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.
Kent
CC BY-SA 4.0
2019-09-17T20:18:50.243Z
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.
gringer
CC BY-SA 4.0
2019-09-17T21:22:01.717Z
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?
Kent
CC BY-SA 4.0
2019-09-17T21:45:53.027Z
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.
gringer
CC BY-SA 4.0
2019-09-16T22:50:48.263Z
2019-09-16T22:50:48.263Z
gringer
<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>

CC BY-SA 4.0
2019-09-17T10:46:43.697Z
2019-09-17T10:46:43.697Z
StupidWolf
<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>

CC BY-SA 4.0
2019-09-17T21:40:23.883Z
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).
gringer
CC BY-SA 4.0
2019-09-17T22:01:27.097Z
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.
Kent
CC BY-SA 4.0
2019-09-17T23:04:59.150Z
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"`.
gringer
CC BY-SA 4.0
2019-09-18T13:16:21.210Z
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.
Kent
<rna-seq><limma><batch-effects><tpm>
TPM or rlog(CPM) for comparing expression?
# 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). 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

**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 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 @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 (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.


---

