Normalization and Scaling

Normalization and Scaling

Single-cell RNA sequencing data is count-based: for each cell, we obtain an integer number of Unique Molecular Identifiers (UMIs) mapped to each gene. These raw counts are not directly comparable between cells because they are strongly influenced by technical factors, most importantly sequencing depth (the total number of molecules captured per cell) and, to a lesser extent, mitochondrial content and cell cycle state. A cell with 10,000 total counts will naturally show higher gene counts than a cell with 2,000 total counts, even if the underlying RNA composition is identical. Normalization is the computational step that corrects for these technical differences so that expression values reflect biology rather than library size artifacts. The most widely used approach, popularized by Scanpy and originally Seurat, is library-size normalization followed by a logarithmic transform. Each cell's counts are scaled by a fixed target (commonly the median count depth or 10,000 counts per cell), and the result is log1p-transformed, i.e., log(1 + x). The logarithm is essential because UMI count distributions are highly skewed; the transform compresses the dynamic range, stabilizes variance, and makes downstream statistical methods such as PCA and differential expression testing behave more reliably. Note that log-transformation is only valid on normalized (non-negative, continuous) values, never on raw integer counts directly. After normalization, scaling (also called z-score standardization) is applied per gene: each gene's expression across all cells is shifted to have a mean of zero and a variance of one. This prevents highly expressed or highly variable genes from dominating principal component analysis purely due to their magnitude. In practice, scaling is typically restricted to highly variable genes (HVGs) and is clipped at a maximum value (commonly 10) to limit the influence of outlier cells on downstream dimensionality reduction. It is important to distinguish normalization from batch correction. Normalization removes technical variation within a dataset arising from sequencing depth, whereas batch correction (e.g., Harmony, scVI, or MNN-based methods) addresses systematic differences between separately processed samples or experiments. Normalization is always performed first; batch correction, if needed, comes later, typically after feature selection and initial dimensionality reduction. The following example demonstrates the standard Scanpy workflow for normalization and scaling: import scanpy as sc # Load a preprocessed AnnData object containing raw counts in adata.X adata = sc.read_h5ad("pbmc3k_raw.h5ad") # 1. Normalize total counts per cell to 10,000, then log1p-transform sc.pp.normalize_total(adata, target_sum=1e4) sc.pp.log1p(adata) # 2. Store the normalized expression for later use (e.g., marker plotting) adata.raw = adata # 3. Identify highly variable genes (HVGs) sc.pp.highly_variable_genes(adata, n_top_genes=2000) # 4. Scale HVGs: zero mean, unit variance, clip extreme values at 10 sc.pp.scale(adata, max_value=10) The second example shows how to implement the same operations manually in NumPy, which is instructive for understanding exactly what these functions do under the hood: import numpy as np # X is a cells x genes matrix of raw UMI counts counts = np.array([ [10, 0, 5, 100], [30, 2, 15, 300], [ 5, 1, 2, 50], ]) # Step 1: library-size normalization to a target of 10,000 target = 1e4 lib_sizes = counts.sum(axis=1, keepdims=True) # per-cell totals normalized = counts / lib_sizes * target # Step 2: log1p transformation log_normalized = np.log1p(normalized) # Step 3: per-gene scaling (z-score) with clipping gene_mean = log_normalized.mean(axis=0) gene_std = log_normalized.std(axis=0) scaled = np.clip((log_normalized - gene_mean) / gene_std, -10, 10) print("Library sizes:", lib_sizes.ravel()) print("Normalized (rounded):", np.round(normalized, 1)) print("Scaled (rounded):", np.round(scaled, 2)) A common pitfall for newcomers is applying these steps in the wrong order or applying them twice. Always confirm whether your input matrix contains raw counts or already-normalized values before calling normalization functions; running normalize_total on log-transformed data produces biologically meaningless results. When in doubt, inspect adata.X for integer values and check whether the object documentation specifies raw counts. Raw UMI counts must be normalized to correct for differences in sequencing depth between cells before any downstream analysis. The standard pipeline is: normalize total counts to a fixed target, apply log1p transformation, select highly variable genes, then scale (z-score) with clipping. Log-transformation stabilizes variance and compresses the skewed count distribution, making PCA and clustering more robust. Scaling ensures each gene contributes comparably to dimensionality reduction, regardless of its average expression level. Normalization is not a substitute for batch correction; the two address different sources of technical variation and are applied at different stages. Keep a copy of normalized, unscaled data (e.g., adata.raw ) so you can visualize meaningful expression values on plots later. Practice exercise: Using the PBMC 3k dataset (available via scv.datasets.pbmc68k_reduced() or the Scanpy tutorial data), write a Python script that loads the raw counts, prints the total count distribution per cell (minimum, median, and maximum), applies sc.pp.normalize_total with target_sum=1e4 followed by sc.pp.log1p , and then compares the expression distribution of the gene LYZ before and after normalization using a histogram. Reflect in a short comment: why does the shape of the distribution change after log-transformation, and what would happen if you mistakenly normalized the data twice?