Why Raw Counts Mislead: Normalization

By now you have aligned your reads and summarized them into a gene-by-sample count matrix with featureCounts or HTSeq. It is tempting to read those integers as expression levels and compare them across genes or across samples directly. Two technical artifacts make that comparison wrong before any biology enters.

Two confounders

The first is library size, or sequencing depth: if one sample yields 27 million assigned reads and another yields 16 million, nearly every gene in the deeper sample carries more counts regardless of biology, so a larger count for the same gene across samples may just reflect depth. The second is gene length: at the same molar abundance a 5 kb transcript is fragmented into more reads than a 1 kb transcript, so within one sample a longer gene looks more expressed than a shorter one even when their true abundances are equal.

Within-sample measures

TPM and FPKM/RPKM divide out gene length as well as a per-sample scaling, so you can compare different genes inside a single sample, and TPM is usually preferred because its per-sample total is fixed at one million while FPKM's is not. CPM divides only by library size with no length term, which makes the same gene roughly comparable across samples for quick QC and clustering. All three are simple per-sample rescalings, and none of them is a valid input to a count-based differential expression test.

library(edgeR)

# gene-by-sample raw counts from featureCounts / HTSeq
counts <- as.matrix(read.delim("counts.tsv", row.names = 1))

# sequencing depth differs across the four samples
colSums(counts)

# CPM = divide each column by its library size, scale to per million
cpm_vals <- cpm(counts)
round(cpm_vals[1:3, ], 1)

   ctrl1    ctrl2   treat1   treat2 
21043918 18755201 27310884 15980433 

          ctrl1  ctrl2 treat1 treat2
ENSG0001   46.0   45.9   45.8   45.9
ENSG0002  304.2  304.6  585.1  576.6
ENSG0003    2.0    1.9    2.0    1.9

Between-sample normalization

For differential expression you need a single scaling factor per sample that survives composition bias, where a handful of very highly expressed genes inflate the total library size and drag simple CPM or TPM off. DESeq2 estimates it with median-of-ratios: take the geometric mean of each gene across all samples, divide each sample's count by its gene's geometric mean, and set the sample's size factor to the median of those ratios, with any gene that has a zero count in any sample dropping out. edgeR's TMM does the analogous job with a trimmed mean of log-ratios against a reference sample, and both rely on the assumption that most genes are not differentially expressed.

library(DESeq2)

coldata <- data.frame(
  condition = factor(c("ctrl", "ctrl", "treat", "treat")),
  row.names = colnames(counts)
)

dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData   = coldata,
                              design    = ~ condition)

# median-of-ratios size factors
dds <- estimateSizeFactors(dds)
sizeFactors(dds)

# counts rescaled by those factors
round(head(counts(dds, normalized = TRUE), 3), 1)

    ctrl1     ctrl2    treat1    treat2 
1.0243117 0.9187446 1.3106842 0.7834529 

          ctrl1   ctrl2  treat1  treat2
ENSG0001  945.0   937.1   955.2   935.6
ENSG0002 6250.1  6218.3 12192.1 11760.8
ENSG0003   42.0    38.1    41.2    38.3