Understanding RNA-Seq Normalization: TPM, FPKM, and Beyond

The Compositional Nature of RNA-Seq Data

A high-throughput RNA-sequencing experiment does not count the absolute number of mRNA molecules inside a living cell. Instead, it measures a random sampling of cDNA fragments competing for finite physical binding space across the lanes of an Illumina flow cell or flow cell nanowells.

Consequently, RNA-Seq count data are fundamentally compositional:

If Gene $A$ increases its expression tenfold, it consumes a larger fraction of total sequencing reads, causing unaffected Genes $B, C, \text{and } D$ to appear downregulated in raw count data simply because fewer total reads remain available to sample them.
Long genes yield more fragmented cDNA pieces than short genes transcribed at the identical cellular abundance.
Samples sequenced with 50 million reads will show proportionally higher raw counts than samples sequenced with 20 million reads, regardless of biological transcription rates.

To derive meaningful biological inferences, bioinformaticians must adjust for these technical confounders through rigorous mathematical normalization.

RPKM and FPKM: The Perils of Non-Constant Denominators

In the early days of short-read transcriptomics (Mortazavi et al., 2008), Reads Per Kilobase of transcript per Million mapped reads (RPKM)—and its paired-end equivalent, Fragments Per Kilobase of transcript per Million mapped fragments (FPKM)—became standard metrics:

$\text{FPKM}_i = \frac{C_i}{\left( \frac{L_i}{10^3} \right) \times \left( \frac{N}{10^6} \right)}$

Where:

$C_i$ is the number of fragments mapped to gene $i$ .
$L_i$ is the exonic length of gene $i$ in base pairs.
$N$ is the total number of mapped fragments in the sequencing library.

The Fatal Flaw of FPKM

While FPKM normalizes for sequencing depth first and gene length second, the sum of normalized FPKM values varies across different sequencing libraries:

$\sum_{i} \text{FPKM}_i \neq \text{Constant}$

Because the denominator differs between biological replicates, you cannot directly compare FPKM values across samples to infer proportional cellular changes. An FPKM of 25 in Sample 1 does not represent the same fraction of cellular transcripts as an FPKM of 25 in Sample 2.

TPM: The Correct Measure of Relative Cellular Abundance

To rectify this mathematical inconsistency, Wagner et al. (2012) proposed Transcripts Per Million (TPM).

TPM reverses the order of operations: it normalizes for gene length first, creating a measure of relative transcript copies, and normalizes for total library abundance second:

$\text{RPK}_i = \frac{C_i}{L_i / 1000}$

$\text{TPM}_i = \frac{\text{RPK}_i}{\sum_j \text{RPK}_j} \times 10^6$

By normalizing by the sum of length-normalized reads ( $\sum_j \text{RPK}_j$ ), the sum of TPM across all genes in every sample is always identically one million:

$\sum_{i} \text{TPM}_i = 1,000,000$

Every 1 TPM unit represents exactly one transcript out of one million cellular transcripts. This makes TPM the correct, reproducible metric for comparing relative expression levels within a sample or across samples for visualization and clustering.

The Library Composition Trap

While TPM is ideal for clustering, heatmaps, and exploratory visualization, neither TPM nor FPKM is valid for statistical differential expression analysis.

Why? Because both assume that total cellular RNA content is invariant across conditions.

Consider a practical biological scenario:

Condition 1: Normal liver tissue.
Condition 2: Liver tissue treated with an oncogenic stimulant that drives massive overexpression of albumin and ribosomal genes, accounting for 60% of all cellular mRNA.

Condition 1 (Baseline)              Condition 2 (Hyper-transcribed)
+-------------------------------+   +-------------------------------+
| Gene A: 20% | Gene B: 20%     |   | Albumin: 60% of all reads    |
| Other Genes: 60%              |   | Gene A: 10% | Gene B: 10%     |
+-------------------------------+   +-------------------------------+
Total mRNA: 1x                      Total mRNA: 3x

In Condition 2, because Albumin consumes 60% of all sequencing reads, Gene $A$ and Gene $B$ will see their relative read proportions cut in half. A standard TPM or FPKM comparison will report that Gene $A$ and Gene $B$ are “significantly downregulated” by 50%, even though their absolute transcription rate inside the cell remained completely unchanged.

DESeq2 Median of Ratios vs EdgeR TMM

To resolve library composition bias, specialized differential expression algorithms (DESeq2 and edgeR) calculate robust scale factors based on invariant housekeeping subsets.

DESeq2: Median of Ratios Method

DESeq2 constructs a pseudo-reference sample by taking the geometric mean of each gene across all libraries:

Calculate Geometric Mean per Gene: $\bar{g}_i = \left( \prod_{j=1}^m K_{ij} \right)^{1/m}$
Calculate Ratio for Gene $i$ in Sample $j$ : $\text{Ratio}_{ij} = \frac{K_{ij}}{\bar{g}_i}$
Determine Sample Size Factor $s_j$ : $s_j = \text{median}_{i} \left( \text{Ratio}_{ij} \right)$

Because the size factor is derived from the median ratio across tens of thousands of genes, it is mathematically immune to extreme outliers (such as a few hyper-expressed oncogenes). Genes that maintain constant basal transcription anchor the scaling factor, restoring true biological fold changes.

Implementing Normalization in Python

Here is a self-contained implementation demonstrating TPM conversion and DESeq2-style size factor estimation from a raw gene count matrix:

import numpy as np
import pandas as pd

def calculate_tpm(counts: pd.DataFrame, gene_lengths_kb: pd.Series) -> pd.DataFrame:
    """
    Computes Transcripts Per Million (TPM) from raw count matrix and gene lengths.
    counts: DataFrame with genes as rows and samples as columns.
    gene_lengths_kb: Series with exonic length in kilobases (bp / 1000).
    """
    # 1. Normalize for gene length (Reads Per Kilobase)
    rpk = counts.divide(gene_lengths_kb, axis=0)
    
    # 2. Normalize for sequencing depth per sample
    scaling_factors = rpk.sum(axis=0) / 1e6
    tpm = rpk.divide(scaling_factors, axis=1)
    
    return tpm

def estimate_deseq2_size_factors(counts: pd.DataFrame) -> pd.Series:
    """
    Computes DESeq2-style median-of-ratios size factors.
    """
    # Filter out genes with zero counts across any sample to compute geometric mean
    non_zero = counts[(counts > 0).all(axis=1)]
    
    # Calculate geometric mean across rows
    log_counts = np.log(non_zero)
    geom_means = np.exp(log_counts.mean(axis=1))
    
    # Compute ratio to geometric mean and take sample median
    ratios = non_zero.divide(geom_means, axis=0)
    size_factors = ratios.median(axis=0)
    
    return size_factors

Decision Matrix: Which Metric to Use?

Analysis Goal	Recommended Metric	Tool / Formula	Why?
Differential Expression	Raw Counts with Size Factors	DESeq2 / edgeR (TMM)	Statistical dispersion models (Negative Binomial) require un-normalized discrete integer counts.
Gene Clustering & Heatmaps	VST or rlog transformed counts	DESeq2 `vst()`	Stabilizes variance across low-count and high-count genes.
Sample-to-Sample Abundance	TPM (Transcripts Per Million)	Salmon / Kallisto	Length-normalized and guarantees equal sample sums.
Cross-condition Comparison	Normalized Counts ( $K_{ij} / s_j$ )	DESeq2 Normalized Counts	Corrects library composition skews and sequencing depth.

Understanding these distinctions prevents researchers from falling into the trap of analyzing FPKM values with linear statistics—a mistake that still accounts for hundreds of irreproducible transcriptomic publications each year.