<p>First off - <strong>do not use FPKM for across-sample normalization</strong>. This is such a common issue that I have a prepared text for it:</p>
<blockquote>
<p>Why not to use FPKM</p>
</blockquote>
<blockquote>
<p>They have a poor sensitivity when looking at low expressed genes [1, 2]. Short transcripts/genes may have falsely over-estimated expression values [5].
They may not normalise adequately across samples [3]. The reason for this is: FPKM are normalised across libraries using a single variable, namely the “total number of reads”. This library size factor estimation can easily be biased by highly expressed, highly variable genes, unlike the method used by the other approaches, which model the library size factor based on summary statistics extracted from all genes. In other words as stated in [4]: “The median ratio method is usually quite robust. The biggest problem for library size correction is not actually the genes with low counts, but the genes with very high counts, which have high variance and disproportionately influence estimators like the total sum. This is why the total sum is a notoriously bad estimator for library size.” Or the same by Lior’s Pachter himself [6].
They have a high False Positive Rate in differential expression analysis [1, 2]. A lot of the genes/transcripts deemed to be differentially expressed are not, mostly due to the effects of point 1), i.e. high artefactual variability in the expression and point 2) i.e. inadequate library size estimation.</p>
</blockquote>
<blockquote>
<p>On the other hand, approaches such as DESeq, DESeq2, vst or voom have adequate between-sample normalisation, have a controlled FDR and can reliably determine differential expression even for lowly expressed genes if there is enough replication or if an adequate log2 fold change cutoff is selected [7]. If you really insist on using cufflinks, consider performing upper quartile normalisation, which corrects for the point 1) issue.</p>
</blockquote>
<blockquote>
<p>References
[1] Dillies M.A., Briefings in Bioinf. 2012
[2] Soneson &amp; Delorenzi, BMC Bioinf. 2013
[3] Lior Pachter’s lecture at CSHL (Cold Spring Harbour Laboratory) from minutes 32 on: <a href="https://www.youtube.com/watch?v=5NiFibnbE8o" rel="nofollow noreferrer">https://www.youtube.com/watch?v=5NiFibnbE8o</a>
[4] From Answer: &quot;Over-correction&quot; in the size-factors of the DESeq2 package:
[5] <a href="http://www.biomedcentral.com/1471-2105/14/370" rel="nofollow noreferrer">http://www.biomedcentral.com/1471-2105/14/370</a>
[6] <a href="https://liorpachter.wordpress.com/2014/04/30/estimating-number-of-transcripts-from-rna-seq-measurements-and-why-i-believe-in-paywall/" rel="nofollow noreferrer">https://liorpachter.wordpress.com/2014/04/30/estimating-number-of-transcripts-from-rna-seq-measurements-and-why-i-believe-in-paywall/</a>
[7] <a href="http://rnajournal.cshlp.org/content/22/6/839" rel="nofollow noreferrer">http://rnajournal.cshlp.org/content/22/6/839</a></p>
</blockquote>
<p>What I recommend, is that you get the raw count data from TCGA, create a count matrix using the DESeq2 package in R, use the function <code>varianceStabilizingTransformation()</code> from DESeq2 to further process your counts and model the networks from there. Your pipeline already mentions <code>voom</code>, which is another good alternative, as is <code>limma</code>. All of these are available on Bioconductor.</p>

CC BY-SA 4.0
2020-10-17T19:41:07.370Z
2020-10-17T19:41:07.370Z
Bastian Schiffthaler
<p>I am trying to collect gene expression data for the point of fitting gene regulatory networks. My background is primarily in computer science and so I am finding the biological literature a bit difficult to penetrate.</p>
<p>Specifically, I want to fit two regulatory networks for several cancers in TCGA from tumor and normal samples. I have queries using the <code>R</code> package <code>TCGABiolinks</code> and was able to follow the tutorial <a href="https://costalab.ukaachen.de/open_data/Bioinformatics_Analysis_in_R_2019/BIAR_D3/handout.html#2_tcga_data" rel="nofollow noreferrer">here</a> up to the normalization part.</p>
<p>So I have three questions;</p>
<ul>
<li>Is there a standard way to normalize gene expression counts? If so, what is it, and is it available as a library in <code>R</code>, <code>python</code> or some other language?</li>
<li>If not, can I just center and scale (z-transform) the samples? I am <strong>not</strong> interested in differential expression, only quantifying <strong>expression</strong> from counts.</li>
</ul>

CC BY-SA 4.0
2020-09-08T08:18:33.433Z
Question is, what do you mean by 'quantifying expression from counts'? TCGA has a standard normalization available, which is FPKM. It normalizes expression by fragment-size and library depth. Some people prefer TPM to FPKM though (it alters the order of operations), so you may want to download the count data and do a TPM normalization using `edgeR`.  For visualization, for TCGA, I personally tend to log2(FPKM) and row-center for heatmaps.
csgroen
CC BY-SA 4.0
2020-09-08T12:10:37.140Z
okay so is the FPKM normalization what is given in the `scaled_estimate` design matrix?
WedgeAntilles
CC BY-SA 4.0
2020-09-08T12:17:29.483Z
Not sure what that matrix is... in TCGAbiolinks or any other way to download from the GDC, you'll have FPKM. It is not scaled though, so probably no.
csgroen
<r>
Standard Way to Preprocess Gene Expression?
# Standard Way to Preprocess Gene Expression?

(Asked by: username_0 on 2020-09-08)

I am trying to collect gene expression data for the point of fitting gene regulatory networks. My background is primarily in computer science and so I am finding the biological literature a bit difficult to penetrate.

Specifically, I want to fit two regulatory networks for several cancers in TCGA from tumor and normal samples. I have queries using the `R` package `TCGABiolinks` and was able to follow the tutorial [here](https://costalab.ukaachen.de/open_data/Bioinformatics_Analysis_in_R_2019/BIAR_D3/handout.html#2_tcga_data) up to the normalization part.

So I have three questions;

* Is there a standard way to normalize gene expression counts? If so, what is it, and is it available as a library in `R`, `python` or some other language?
* If not, can I just center and scale (z-transform) the samples? I am **not** interested in differential expression, only quantifying **expression** from counts.

**Tags:** r
> **Comments:**
> **username_1**: Question is, what do you mean by 'quantifying expression from counts'? TCGA has a standard normalization available, which is FPKM. It normalizes expression by fragment-size and library depth. Some people prefer TPM to FPKM though (it alters the order of operations), so you may want to download the count data and do a TPM normalization using `edgeR`.  For visualization, for TCGA, I personally tend to log2(FPKM) and row-center for heatmaps.  
> **username_0**: okay so is the FPKM normalization what is given in the `scaled_estimate` design matrix?  
> **username_1**: Not sure what that matrix is... in TCGAbiolinks or any other way to download from the GDC, you'll have FPKM. It is not scaled though, so probably no.  

---

## 1 answers

### Answer by username_2 on 2020-10-17. Score: 2

First off - **do not use FPKM for across-sample normalization**. This is such a common issue that I have a prepared text for it:

> Why not to use FPKM

> They have a poor sensitivity when looking at low expressed genes [1, 2]. Short transcripts/genes may have falsely over-estimated expression values [5].
> They may not normalise adequately across samples [3]. The reason for this is: FPKM are normalised across libraries using a single variable, namely the “total number of reads”. This library size factor estimation can easily be biased by highly expressed, highly variable genes, unlike the method used by the other approaches, which model the library size factor based on summary statistics extracted from all genes. In other words as stated in [4]: “The median ratio method is usually quite robust. The biggest problem for library size correction is not actually the genes with low counts, but the genes with very high counts, which have high variance and disproportionately influence estimators like the total sum. This is why the total sum is a notoriously bad estimator for library size.” Or the same by Lior’s Pachter himself [6].
> They have a high False Positive Rate in differential expression analysis [1, 2]. A lot of the genes/transcripts deemed to be differentially expressed are not, mostly due to the effects of point 1), i.e. high artefactual variability in the expression and point 2) i.e. inadequate library size estimation.

> On the other hand, approaches such as DESeq, DESeq2, vst or voom have adequate between-sample normalisation, have a controlled FDR and can reliably determine differential expression even for lowly expressed genes if there is enough replication or if an adequate log2 fold change cutoff is selected [7]. If you really insist on using cufflinks, consider performing upper quartile normalisation, which corrects for the point 1) issue.

> References
> [1] Dillies M.A., Briefings in Bioinf. 2012
> [2] Soneson & Delorenzi, BMC Bioinf. 2013
> [3] Lior Pachter’s lecture at CSHL (Cold Spring Harbour Laboratory) from minutes 32 on: <https://www.youtube.com/watch?v=5NiFibnbE8o>
> [4] From Answer: "Over-correction" in the size-factors of the DESeq2 package:
> [5] <http://www.biomedcentral.com/1471-2105/14/370>
> [6] <https://liorpachter.wordpress.com/2014/04/30/estimating-number-of-transcripts-from-rna-seq-measurements-and-why-i-believe-in-paywall/>
> [7] <http://rnajournal.cshlp.org/content/22/6/839>

What I recommend, is that you get the raw count data from TCGA, create a count matrix using the DESeq2 package in R, use the function `varianceStabilizingTransformation()` from DESeq2 to further process your counts and model the networks from there. Your pipeline already mentions `voom`, which is another good alternative, as is `limma`. All of these are available on Bioconductor.


---

