Normalization

Two samples in a count matrix were almost never sequenced to the same depth. The deeper library has more reads at nearly every gene, and because that difference belongs to the run rather than the biology, the raw counts are not directly comparable between samples. DESeq2 corrects for depth with size factors, one number per sample, dividing each sample’s counts by its factor so that a count reflects expression rather than the depth of its library.

The size factor

The reference is the geometric mean of each gene’s counts across samples, one value per gene and the same for every sample. A sample’s size factor is the median, over genes, of its count divided by that reference. If one sample was sequenced twice as deep as the rest, every ratio for that sample doubles, and its size factor comes out near 2.

library(DESeq2)
library(airway)

data(airway)
dds <- DESeqDataSet(airway, design = ~ cell + dex)
dds$dex <- relevel(dds$dex, ref = "untrt")

# Estimate size factors.
dds <- estimateSizeFactors(dds)

# One size factor per sample.
sizeFactors(dds)

SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520
 1.0236476  0.8961667  1.1794861  0.6700538  1.1776714  1.3990365  0.9207787
SRR1039521
 0.9445141

All eight factors here sit near 1, so the libraries were sequenced to a similar depth, with SRR1039513 the shallowest at 0.6700538 and SRR1039517 the deepest at 1.3990365. Dividing each sample’s counts by its own factor brings the samples onto a common scale.

Raw and normalized counts

The first gene in the matrix has 679 raw reads in SRR1039508 and 448 in SRR1039509, and part of that difference is depth rather than expression. Division by the size factors gives 663.3142 and 499.9070, and on that scale the two samples are directly comparable.

# Raw counts for one gene, across samples.
head(counts(dds)[1, ])

# Normalized counts for the same gene.
head(counts(dds, normalized = TRUE)[1, ])

SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517
       679        448        873        408       1138       1047
SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517
  663.3142   499.9070   740.1528   608.9063   966.3137   748.3722

The differential expression test runs on the raw counts, with the size factors inside the model as offsets, because the negative binomial model needs integer counts. The normalized matrix is for plotting and for sanity checks.

The variance-stabilizing transform

Depth is not the only problem with plotting counts. Variance grows with the mean in count data, so on the count scale the genes with the highest means dominate a sample distance or a clustering. The variance-stabilizing transform removes that dependence, and on the transformed matrix a gene’s contribution to a distance no longer grows with its count.

# Variance-stabilizing transform, using the size factors.
vsd <- vst(dds, blind = FALSE)

# The transformed matrix.
dim(assay(vsd))
head(assay(vsd)[, 1:3], 3)

[1] 63677     8
                SRR1039508 SRR1039509 SRR1039512
ENSG00000000003   9.742074   9.430420   9.867627
ENSG00000000005   6.681995   6.681995   6.681995
ENSG00000000419   9.333669   9.581707   9.486145

With blind = FALSE the dispersion trend is estimated with the design in the model, so the transform does not mistake the condition effect for noise. Once the QC step is done, that is the right setting. The QC page used blind = TRUE for an unsupervised look at the samples. Either setting returns a matrix on the same scale; all that changes is what informed the estimate.

What normalization does not do

A size factor corrects for depth and for nothing else. Gene length stays in the counts, so the normalized matrix still cannot compare one gene against another within a sample the way a TPM can. For a comparison between conditions the omission does not matter, since length is a constant of the gene and cancels out of the fold change, and DESeq2 leaves it unadjusted. The differential expression page picks up the size-factor-adjusted object and fits the model.