bio-expression-matrix-normalization

What this skill does

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

| Biology | Mechanism | Symptom |
|---------|-----------|---------|
| MYC amplification | MYC drives global 2-3x transcriptional amplification | Scaling factors absorb the global shift; reported LFCs muted |
| Apoptosis / cell death | Massive transcriptional shutdown | Surviving (often mitochondrial) transcripts appear up-regulated |
| Viral host shutoff (HSV-1, vaccinia) | Host mRNA degraded by viral nucleases | Apparent up-regulation of non-degraded host genes |
| Transcription / splicing inhibitors (DRB, flavopiridol) | Pol II elongation blocked | Similar global shutdown signature |
| Cell-cycle synchronization | G1 vs S vs M differ in total RNA per cell | Per-cell RNA shifts; library-size-only normalization mis-corrects |
| Prokaryotic stress | Bacteria rewire large fractions of transcriptome | >50% of genes truly DE; trimmed mean is no longer "background" |

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

| Task | Method | Tool | Input |
|------|--------|------|-------|
| DE testing (DESeq2) | RLE / median-of-ratios | `DESeq()` internally | Raw integer counts |
| DE testing (edgeR) | TMM / TMMwsp | `normLibSizes()` internally | Raw integer counts |
| DE testing (limma-voom) | TMM + voom precision weights | `normLibSizes()` then `voom()` | Raw integer counts |
| PCA, heatmaps, clustering (n > 30) | VST | `vst(dds, blind = FALSE)` | DESeqDataSet |
| PCA, heatmaps, clustering (n < 30, library sizes vary >4x) | rlog | `rlog(dds, blind = FALSE)` | DESeqDataSet |
| PCA, heatmaps (edgeR/limma) | log-CPM | `cpm(y, log=TRUE, prior.count=2)` | DGEList |
| WGCNA | VST or log-CPM | `vst(blind=FALSE)` | DESeqDataSet |
| GSVA / ssGSEA | log2(TPM+1) or VST | precomputed | TPM or DESeqDataSet |
| ML biomarker model | VST | `vst(blind=FALSE)` | DESeqDataSet |
| Cross-sample expression reporting | DESeq2 normalized counts | `counts(dds, normalized=TRUE)` | DESeqDataSet |
| Within-sample gene ranking | TPM | quantification tool output | -- |
| Single-cell normalization | scran deconvolution | `computeSumFactors` + `logNormCounts` | SingleCellExperiment |
| Single-cell, simple/exploratory | log1p(CPM-like) | `sc.pp.normalize_total + sc.pp.log1p` | AnnData |
| Cross-sample but composition-shifted (majority-DE biology) | Spike-in (ERCC) or `controlGenes=` | DESeq2/edgeR with offset | Raw counts + controls |
| GC-content bias (cross-platform integration) | cqn or EDASeq | `cqn()` returns offset | Counts + per-gene GC + length |

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

```r
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:

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

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

| Scenario | Use |
|----------|-----|
| Standard bulk RNA-seq | Default `type='ratio'` |
| Zero-heavy (single-cell, sparse) | `type='poscounts'` |
| Very small libraries | `type='iterate'` |
| Known stable reference genes | `controlGenes=stable_idx` |
| Majority-DE biology (stress, MYC, viral) | `controlGenes=` or spike-in SBN |

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

```r
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:

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

```r
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 panel