Chapter 4: Normalization, Feature Selection & Variance Stabilization

Chapter 4: Normalization, Feature Selection & Variance Stabilization

Single-cell RNA sequencing suffers from extreme technical variability: sequencing depth differences, capture efficiency variability between droplets, and dropout effects (false zeros where an expressed mRNA was not captured). Normalization decouples technical sequencing depth from true biological gene expression.

1. Classical Global Scaling (LogNormalize)

For each gene \(i\) in cell \(j\):

y_{ij} = \ln\left( rac{x_{ij}}{\sum_{k} x_{kj}} imes 10,000 + 1 ight)

# Standard Seurat LogNormalize workflow
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)

# Identify 2,000 Highly Variable Genes (vst method)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

# Scale data and regress out unwanted nuisance sources (e.g. percent.mt)
pbmc <- ScaleData(pbmc, vars.to.regress = c("percent.mt"))

2. Regularized Negative Binomial Regression (SCTransform v2)

Developed by Hafemeister & Satija, SCTransform models single-cell counts as a generalized linear model with negative binomial error distribution. It calculates Pearson residuals, effectively removing sequencing depth while preserving biological heterogeneity without manual scaling or pseudocounts.

# Running SCTransform v2 in Seurat v5
pbmc <- SCTransform(
    pbmc, 
    vst.flavor = "v2", 
    vars.to.regress = "percent.mt", 
    variable.features.n = 3000, 
    verbose = FALSE
)

# Active assay is now automatically set to "SCT"
# SCT assay contains Pearson residuals in the 'scale.data' slot for PCA/UMAP
DefaultAssay(pbmc) <- "SCT"

3. Highly Variable Gene (HVG) Selection

Most of the 30,000+ genes in the genome are stably expressed housekeeping genes or low-abundance noise. Performing clustering across all genes dilutes the biological signal with technical noise. HVG selection fits a polynomial curve modeling the mean-variance relationship and picks the top 2,000–3,000 genes with excess biological variance.

# Scanpy HVG identification with sample batch awareness
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

# Compute highly variable genes using Seurat or cell_ranger flavor
sc.pp.highly_variable_genes(
    adata, 
    n_top_genes=2500, 
    subset=False, 
    flavor='seurat_v3', 
    batch_key='batch' # Prevent picking genes variable only due to batch differences
)

# Plot dispersion vs mean expression
sc.pl.highly_variable_genes(adata)

4. Cell Cycle Scoring & Regression

In dividing cells, proliferation genes (cyclins, CDKs, DNA replication machinery) dominate principal components, clustering cells by cell cycle phase (G1, S, G2/M) rather than biological cell identity.

# Scoring cell cycle phases in Seurat
s.genes <- cc.genes$s.genes
g2m.genes <- cc.genes$g2m.genes

pbmc <- CellCycleScoring(pbmc, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)

# If cell cycle is a technical confounder, regress out the difference
pbmc$CC.Difference <- pbmc$S.Score - pbmc$G2M.Score
pbmc <- SCTransform(pbmc, vars.to.regress = c("percent.mt", "CC.Difference"))