Bio Expression Matrix Normalization

Version Compatibility

Reference examples tested with: DESeq2 1.42+, edgeR 4.0+, limma 3.58+, pandas 2.2+, numpy 1.26+, scanpy 1.10+, scran 1.30+, scater 1.30+, EDASeq 2.36+, cqn 1.48+

Before using code patterns, verify installed versions match. If versions differ:

• R: packageVersion('<pkg>') then ?function_name to verify parameters
• Python: pip show <package> then help(module.function) to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Expression Matrix Normalization

"Normalize my counts for X" -> Pick a method that matches the downstream task (raw counts for DE; TMM/RLE-scaled CPM for cross-sample comparison; VST/rlog for visualization/ML; scran for single-cell) and the data structure (bulk vs single-cell, zero-heavy or not, length-biased or not).

The Single Most Important Modern Insight -- TMM/RLE assume most genes are NOT DE; that assumption is wrong in MYC, apoptosis, viral host shutoff, and prokaryotic stress

TMM (Robinson & Oshlack 2010 Genome Biol 11:R25) and RLE/median-of-ratios (Anders & Huber 2010 Genome Biol 11:R106) both rely on the assumption that the MAJORITY of genes are unchanged between samples. The scaling factor is computed from a trimmed reference. When that assumption holds (most bulk RNA-seq), both methods are excellent.

When the assumption fails catastrophically:

Detection: MA plot shows the bulk cloud clearly shifted off zero; reported fold changes don't match qPCR or Western for known-DE genes.

Fix: ERCC spike-in normalization (Jiang 2011 Genome Res 21:1543) -- 96 synthetic RNAs spiked proportional to cell number; use spike-in counts as the size factor. Or: controlGenes= with a curated stable housekeeping set. The "most genes not DE" assumption is checkable; check it on every dataset that might violate it.

A second insight: lengthScaledTPM from tximport is NOT TPM. It is a count-scale matrix (sums to library size) with length-bias removed, designed as input to limma-voom (which cannot accept offsets). The "TPM" in the option name has misled many users into reporting these values as normalized abundance. See expression-matrix/counts-ingest for the full countsFromAbundance decision tree.

A third insight: VST/rlog values are NEVER input to DE. They are for PCA, heatmaps, clustering, ML features. DE tools model the count distribution directly; passing VST/rlog values silently violates their assumptions.

Normalization Decision Table

Between-Sample Normalization

RLE / Median of Ratios (DESeq2)

Goal: Estimate per-sample size factors that correct for library size and composition bias.

Approach: Geometric mean per gene across samples forms a pseudo-reference; per-sample size factor is the median ratio of (sample count / reference count). Median is robust to a minority of DE genes.

library(DESeq2)

dds <- DESeqDataSetFromMatrix(counts, coldata, design = ~ condition)
dds <- estimateSizeFactors(dds)
sizeFactors(dds)

norm_counts <- counts(dds, normalized = TRUE)

Size factor interpretation: 1.2 means the sample has 20% more sequencing depth (after composition adjustment) than the reference.

The geometric mean is undefined for any gene with a zero in any sample, so RLE type='ratio' (default) silently EXCLUDES genes with any zero from the reference. For single-cell or zero-heavy data, use:

dds <- estimateSizeFactors(dds, type = 'poscounts')

poscounts uses only positive entries per gene, salvaging the geometric mean.

TMM / TMMwsp (edgeR)

Goal: Compute normalization factors that account for composition bias via trimmed mean of M-values.

Approach: Select reference; compute gene-wise log-ratios (M-values) and average expression (A-values); trim extremes; weighted mean scaling factor.

library(edgeR)

y <- DGEList(counts = counts, group = coldata$condition)
y <- normLibSizes(y)
y$samples$norm.factors

normLibSizes() is the v4 canonical name (was calcNormFactors() in v3; same function). method='TMM' remains the documented default; method='TMMwsp' (TMM with singleton pairing) is an alternative for samples with many zeros and is the preferred choice for sparse / single-cell-pseudobulk data. Pass method= explicitly for reproducibility. Both old name and method still work.

TMM defaults: logratioTrim = 0.3 (trim top 30% and bottom 30% of M values), sumTrim = 0.05. The factor enters the GLM as part of the offset; counts are NOT divided.

For visualization:

log_cpm <- cpm(y, log = TRUE, prior.count = 2)
cpm_vis <- cpm(y, normalized.lib.sizes = TRUE)

prior.count = 2 is the modern edgeR default for log=TRUE. Smaller priors (0.25) leave low-count log values noisy; larger priors (5-10) shrink them further toward zero. limma-trend assumes the log-CPM uses cpm log=TRUE.

Upper Quartile

y <- normLibSizes(y, method = 'upperquartile')

Per-sample 75th percentile of non-zero counts (Bullard et al. 2010 BMC Bioinformatics 11:94). Robust to the "most genes unchanged" assumption -- useful when TMM/RLE fail. Less common in modern practice; usable for microRNA-seq or targeted panels where TMM doesn't have enough genes to trim.

Spike-In Normalization (ERCC)

ercc_idx <- grep('^ERCC-', rownames(counts))
y <- DGEList(counts = counts[-ercc_idx, ], group = coldata$condition)
ercc_lib <- colSums(counts[ercc_idx, ])
y$samples$norm.factors <- ercc_lib / mean(ercc_lib)

Spike-in based normalization (SBN) uses ERCC counts as the per-sample size factor proxy, decoupling normalization from endogenous transcript composition. Required for MYC, viral, or other majority-DE biology. Caveat: spike-ins have technical variance; the spike-in:endogenous ratio depends on input cell number being accurately measured (often violated).

Within-Sample Normalization (TPM, FPKM)

import pandas as pd

def counts_to_tpm(counts, gene_lengths):
    '''Convert raw counts to TPM. gene_lengths in bp.'''
    rate = counts.div(gene_lengths / 1000, axis=0)
    tpm = rate.div(rate.sum(axis=0), axis=1) * 1e6
    return tpm

counts_to_tpm <- function(counts, gene_lengths) {
    rate <- counts / (gene_lengths / 1000)
    t(t(rate) / colSums(rate)) * 1e6
}

TPM's sum-to-1e6-per-sample makes cross-sample comparisons dimensionally coherent but does NOT fix composition shifts -- TPM remains a COMPOSITIONAL value. Under massive global changes, TPM destroys the absolute-scale signal.

Wagner GP, Kin K, Lynch VJ 2012 Theory Biosci 131:281 is the canonical "why TPM > RPKM" reference. Their headline: "average FPKM varies between samples even for the same genome."

CRITICAL: do NOT use TPM as input to DESeq2 / edgeR / limma-voom. They model the count distribution directly. TPM is for: within-sample ranking, gene-set comparison within sample, deconvolution (CIBERSORT, EPIC, MCP-counter all require TPM), and reporting expression levels.

Variance-Stabilizing Transformations (VST, rlog)

Goal: Transform counts so variance is approximately constant across the mean, suitable for PCA, heatmaps, clustering, ML features.

Approach: Fit dispersion-mean trend and apply variance-stabilizing function.

library(DESeq2)

vsd <- vst(dds, blind = FALSE)
rld <- rlog(dds, blind = FALSE)

vst_matrix <- assay(vsd)

blind=TRUE (vst default) re-estimates dispersions IGNORING the design -- appropriate for unbiased QC ("are samples consistent independent of design?"). blind=FALSE uses fitted dispersions -- appropriate for downstream visualization after the design is settled. Modern DESeq2 vignette recommends blind=FALSE for any plot AFTER the model is fit.

vst() uses 1000 most-variable genes by default to fit the dispersion trend. With <1000 genes after filtering, set nsub lower.

CRITICAL: VST/rlog values are for visualization only. NEVER pass them to DESeq() or glmQLFit(). DE tools require raw counts.

log-CPM (edgeR / limma)

log_cpm <- cpm(y, log = TRUE, prior.count = 2)

import numpy as np
def log_cpm(counts, prior_count=2):
    lib_sizes = counts.sum(axis=0)
    cpm_vals = (counts + prior_count) / (lib_sizes + 2 * prior_count) * 1e6
    return np.log2(cpm_vals)

prior.count = 2 (edgeR modern default) shrinks low-count log values toward zero, reducing visual artifacts from low-expression noise. For statistical use (limma-trend), this default is appropriate; for heatmaps and PCA, larger priors (3-5) further dampen noise.

voom uses 0.5 added per cell with (lib.size + 1) denominator -- not strictly equivalent to cpm(log=TRUE, prior.count=0.5). Its mean-variance trend handles low-count noise via per-observation weights. Do not confuse voom's internal pre-log handling with cpm(log=TRUE).

GC Content and Length Bias (cqn, EDASeq)

Standard TMM/RLE do NOT correct sample-specific GC content bias or gene length bias. These biases arise from library preparation (fragmentation, PCR) and create systematic differences in read coverage correlated with gene properties.

Most affected: GSEA. Sample-specific gene-length bias causes recurrent false positives.

EDASeq (Risso, Schwartz, Sherlock, Dudoit 2011 BMC Bioinformatics 12:480) applies within-lane GC normalization then between-lane normalization:

library(EDASeq)

fd <- data.frame(gc = gene_gc, length = gene_lengths, row.names = rownames(counts))
data <- newSeqExpressionSet(as.matrix(counts), featureData = fd, phenoData = coldata)

data_norm <- withinLaneNormalization(data, 'gc', which = 'full')
data_norm <- betweenLaneNormalization(data_norm, which = 'full')
normalized_counts <- counts(data_norm)

cqn (Hansen, Irizarry, Wu 2012 Biostatistics 13:204) fits a smooth conditional-quantile model on log2 expression as a function of GC and length; the offset enters the GLM directly:

library(cqn)
cqn_res <- cqn(counts, x = gene_gc, lengths = gene_lengths)
y$offset <- cqn_res$glm.offset

Use when comparing samples sequenced on different platforms or with different library prep chemistries. Run BEFORE TMM/RLE (or supply the cqn offset to edgeR/DESeq2 directly).

Single-Cell Normalization

scran Deconvolution

Goal: Per-cell size factors that handle the high zero-inflation of single-cell data without violating TMM/RLE assumptions.

Approach: Pool cells in overlapping windows; compute pool-level size factors; deconvolve back to individual cells. Pre-cluster to avoid mixing cell types with very different transcriptome sizes.

library(scran)
library(scater)

clusters <- quickCluster(sce)
sce <- computeSumFactors(sce, clusters = clusters)
sce <- logNormCounts(sce)

quickCluster provides a rough partition; without it, mixing a cell with 200 detected genes (resting T cell) and 5000 detected genes (activated macrophage) produces wrong size factors.