Chapter 4: Statistical Foundations of Sequencing Counts (Overdispersion)
Read Counts are Random Variables
If you take a single tube of RNA, prepare two biological libraries, and sequence them under identical conditions, you will not get the same count for any gene. A gene that yields 450 reads in one sample might yield 420 in another. This is not measurement error; it is biological variation and random sampling. A statistical model is not a mathematical convenience—it is an explicit statement about where variation originates.
Core Concepts
In genomics, primary data arrives as discrete counts—integer values representing the number of sequencing reads mapped to a gene or genomic interval. A gene is observed \(0, 1, 2,\) or \(500\) times; it cannot be observed \(450.7\) times.
Why can you not apply standard Gaussian (normal distribution) statistics—such as Student’s \(t\)-test or classic ANOVA—directly to raw, un-transformed sequencing counts?
The core statistical objection is not merely that counts cannot be negative (a normal distribution can approximate positive counts well when far from zero). The fundamental objection is mean-variance heteroscedasticity: in sequencing count data, variance scales non-linearly with the mean. Highly expressed genes exhibit vastly larger absolute variance than lowly expressed genes, violating the normal distribution’s requirement of homoscedasticity (equal variance across comparison groups).
Gaussian-based statistical frameworks (such as limma-voom or DESeq2’s Variance Stabilizing Transformation [VST]) can be applied successfully after transforming counts into continuous space and modeling the mean-variance relationship explicitly. But applying Gaussian models directly to un-transformed, raw counts leads to severely biased inference.
The Poisson Model: Random Sampling Baseline
Sequencing is fundamentally a random sampling process. A flow cell extracts a tiny subset of millions of cDNA molecules from a library containing billions. When drawing rare items independently at a constant rate, the counts are described by the Poisson distribution:
\[P(Y = y) = \frac{\lambda^y e^{-\lambda}}{y!}\]
where \(\lambda\) represents the expected rate (mean count). The defining constraint of the Poisson distribution is that its variance equals its mean:
\[\text{E}[Y] = \lambda \quad \text{and} \quad \text{Var}(Y) = \lambda\]
I suggest you understand what the Poisson model represents: Poisson provides a baseline theoretical model for random sampling noise.
When you sequence technical replicates—for example, taking the exact same tube of RNA library and splitting it across three sequencing lanes—the Poisson model fits the observed count variation very well. The only source of variation is the random chance of which cDNA fragment lands on the flow cell surface.
Technical vs. Biological Replicates
I recommend you memorize this fundamental distinction:
Technical replicates tell you how stable the measurement instrument is.
Biological replicates tell you how variable the underlying biological population is.
Differential expression inference requires biological replicates (e.g., sequencing liver tissue from 3 distinct mice or tumor biopsies from 3 human patients). In biological replicates, true individual variation is layered on top of technical sampling noise.
Library Size Factors: The Compositional Equation
Before analyzing count variation across samples, you must account for a fundamental experimental variable: sequencing depth (library size).
Suppose Sample 1 was sequenced to a depth of 10 million reads, while Sample 2 was sequenced to a depth of 30 million reads. A raw count of 1,000 in Sample 2 does not indicate higher biological expression than a raw count of 500 in Sample 1.
In quantitative omics, the observed count is governed by a compositional equation:
\[\text{Observed Count} = (\text{Biological Abundance}) \times (\text{Library Size Factor}) + \text{Sampling Noise} + \text{Biological Variability}\]
Software like DESeq2 computes a relative size factor (\(s_j\)) for each sample \(j\) (typically using the median ratio method across all genes). Dividing raw counts by size factors yields normalized counts, ensuring that variation across samples reflects true biological differences rather than total flow cell loading.
Overdispersion and the Negative Binomial Model
Consider a real-world transcriptomic measurement across three biological replicates of healthy mouse liver tissue normalized for library size:
Mouse 1: \(150\) reads
Mouse 2: \(300\) reads
Mouse 3: \(450\) reads
Computing sample statistics for this gene:
\[\text{Sample Mean } (\hat{\mu}) = \frac{150 + 300 + 450}{3} = 300\]
\[\text{Sample Variance } (s^2) = \frac{(150-300)^2 + (300-300)^2 + (450-300)^2}{3 - 1} = \frac{22,500 + 0 + 22,500}{2} = 22,500\]
Under a Poisson model (\(\text{Var} = \mu\)), a gene with a mean of \(300\) should have a variance of \(300\). In reality, the observed sample variance (\(22,500\)) is 75 times larger than the mean.
This phenomenon is called overdispersion (\(\text{Var}(Y) > \text{E}[Y]\)). It occurs because individual organisms differ in genetic background, environmental history, cell-type composition, and transcriptional bursting.
Intuitive Definition of Dispersion (\(\alpha\))
I suggest you think of mean and dispersion intuitively:
The mean (\(\mu\)) tells you where the counts are centered.
The dispersion (\(\alpha\)) tells you how much counts can wander around that mean.
In the Poisson model, wandering is strictly constrained by the mean (\(\text{Var} = \mu\)).
In the Negative Binomial (NB) distribution, biological variability gives counts extra room to wander:
\[\text{Var}(Y) = \mu + \alpha \mu^2\]
When \(\alpha = 0\): The quadratic term vanishes, and the model reverts to the Poisson baseline (\(\text{Var} = \mu\)).
When \(\alpha > 0\): The quadratic term \(\alpha \mu^2\) explicitly models biological overdispersion across individual organisms.
Figure Placeholder: [Figure 4.1: Poisson vs. Negative Binomial Probability Distributions. Panel A shows Poisson distributions at low (\(\lambda=10\)) and high (\(\lambda=100\)) means, demonstrating how variance strictly equals the mean. Panel B overlays Negative Binomial distributions with increasing dispersion parameters (\(\alpha=0.01, 0.1, 0.5\)), showing how overdispersion dramatically widens the distribution tails while holding the mean constant at \(\mu=100\).]
To model both technical sampling noise and biological variation, we model the true biological rate \(\lambda\) across replicates not as a fixed constant, but as a random variable following a Gamma distribution:
\[\lambda \sim \text{Gamma}\left(k = \frac{1}{\alpha}, \theta = \alpha \mu\right)\]
Conditional on \(\lambda\), the sequencing read count \(Y\) follows a Poisson distribution:
\[Y \mid \lambda \sim \text{Poisson}(\lambda)\]
Integrating out the unobserved biological rate \(\lambda\) yields the marginal distribution of \(Y\):
\[P(Y = y) = \int_{0}^{\infty} P(Y = y \mid \lambda) f(\lambda) \, d\lambda = \int_{0}^{\infty} \frac{\lambda^y e^{-\lambda}}{y!} \cdot \frac{\lambda^{k-1} e^{-\lambda/\theta}}{\theta^k \Gamma(k)} \, d\lambda\]
Evaluating this integral produces the probability mass function of the Negative Binomial distribution:
\[P(Y = y) = \frac{\Gamma(y + 1/\alpha)}{y! \, \Gamma(1/\alpha)} \left( \frac{1/\alpha}{1/\alpha + \mu} \right)^{1/\alpha} \left( \frac{\mu}{1/\alpha + \mu} \right)^y\]
Using the law of total variance:
\[\text{Var}(Y) = \text{E}_\lambda[\text{Var}(Y \mid \lambda)] + \text{Var}_\lambda(\text{E}[Y \mid \lambda]) = \text{E}_\lambda[\lambda] + \text{Var}_\lambda(\lambda) = \mu + \alpha \mu^2\]
Conceptual Demonstration: Overdispersion vs. Statistical Evidence
Let us examine a conceptual example demonstrating how neglecting overdispersion distorts statistical uncertainty.
Suppose we compare 3 Control mice to 3 Treated mice for Gene X (normalized for library size):
Control counts: \(150, 300, 450 \implies \hat{\mu}_{\text{ctrl}} = 300, \, s^2_{\text{ctrl}} = 22,500\)
Treated counts: \(600, 900, 1200 \implies \hat{\mu}_{\text{trt}} = 900, \, s^2_{\text{trt}} = 90,000\)
Fold Change: \(\frac{900}{300} = 3.0\) (a 3-fold upregulation).
Naive Poisson Estimation
If a naive algorithm assumes a Poisson model (\(\text{Var} = \mu\)), it calculates standard errors assuming variance equals mean (\(\text{Var}_{\text{ctrl}}=300, \text{Var}_{\text{trt}}=900\)):
\[\text{SE}_{\text{Poisson}} = \sqrt{\frac{300}{3} + \frac{900}{3}} = \sqrt{100 + 300} = 20\]
\[\text{Z-score}_{\text{Poisson}} = \frac{900 - 300}{20} = 30.0 \implies p < 10^{-190}\]
The naive Poisson model reports a near-zero \(p\)-value (\(p < 10^{-190}\)). A naive researcher might declare Gene X a massive biological discovery. But a tiny \(p\)-value does not rescue a flawed statistical model.
Negative Binomial Estimation
Incorporating true biological dispersion (\(\alpha \approx 0.2467\)):
\[\text{Var}_{\text{ctrl}} = 300 + 0.2467(300)^2 = 22,500 \quad \implies \quad \text{Var}_{\text{trt}} = 900 + 0.2467(900)^2 = 200,727\]
\[\text{SE}_{\text{NB}} = \sqrt{\frac{22,500}{3} + \frac{200,727}{3}} \approx \sqrt{7,500 + 66,909} = \sqrt{74,409} \approx 272.78\]
\[\text{Z-score}_{\text{NB}} = \frac{900 - 300}{272.78} = 2.20 \implies p \approx 0.028\]
Accounting for biological overdispersion shifts the \(p\)-value from \(10^{-190}\) to \(0.028\).
I ask you to reflect on the core biological lesson: Effect size (fold change) and statistical evidence are not the same thing. A 3-fold expression change can represent weak statistical evidence if biological variability across individuals is large (\(\alpha = 0.25\)). Conversely, a modest 1.5-fold change can provide strong statistical evidence if counts across individuals are tightly constrained.
Dispersion Estimation in Practice: Empirical Bayes Shrinkage
In typical transcriptomic experiments, biological sample sizes are small (\(N = 3\) to \(5\) per group). Estimating an independent dispersion parameter \(\alpha_g\) for each of 20,000 genes from only 3 samples yields noisy, unreliable gene-wise estimates (\(\hat{\alpha}_g^{\text{MLE}}\)).
Modern packages (DESeq2, edgeR) solve small-\(N\) instability using Empirical Bayes MAP Shrinkage:
Figure Placeholder: [Figure 4.2: Empirical Bayes Dispersion Shrinkage in DESeq2. A scatter plot showing mean log normalized counts on the x-axis versus dispersion estimates on the y-axis for ~20,000 genes. Black dots represent raw gene-wise Maximum Likelihood Estimates (MLEs). A red fitted line shows the global mean-dispersion trend (\(\alpha(\mu) = a_1/\mu + a_0\)). Blue arrows show gene-wise estimates being shrunk toward the red trend line to yield final Maximum A Posteriori (MAP) estimates (blue circles). High-dispersion outliers whose raw values sit far above the trend line are not shrunk, preventing under-estimation of variance.]
Gene-wise Maximum Likelihood Estimation (MLE): Compute a raw dispersion estimate \(\hat{\alpha}_g^{\text{MLE}}\) for each gene \(g\).
Mean-Dispersion Trend Fitting: Fit a smooth curve relating mean count \(\mu\) to dispersion \(\alpha\) across all 20,000 genes: \[\alpha(\mu) = \frac{a_1}{\mu} + a_0\] With few samples per gene, individual dispersion estimates are noisy; fitting a trend allows the model to borrow information across genes sharing similar mean expression levels.
Empirical Bayes Maximum A Posteriori (MAP) Shrinkage: Shrink individual gene-wise estimates toward the fitted trend line. Genes with typical variance are shrunk strongly toward the curve, stabilizing variance estimates for differential testing.
- Outlier Handling: Genes with raw dispersion estimates that vastly exceed the trend line (e.g., highly variable surface receptors or immune genes) are flagged as dispersion outliers. Their dispersion is not shrunk downward, ensuring that statistical testing does not under-estimate variance and produce false-positive calls.
Biological Interpretation
Gene Variability and Statistical Power
Gene dispersion is not a fixed universal constant; it varies across biological processes and experimental designs.
Genes mediating foundational cellular survival (such as structural proteins or core metabolic enzymes) often exhibit tight expression control across healthy individuals. Conversely, genes responding to environmental stimuli, immune signaling, or stress exhibit substantial inter-individual variability.
I suggest you understand the practical consequence: Genes with high biological dispersion require larger sample sizes to achieve equivalent statistical power.
There is no universal “magic sample size” (\(N=3\) vs \(N=10\)). Required biological sample size depends on four factors:
Expected biological effect size (fold change).
Gene dispersion (\(\alpha\)).
Sequencing read depth.
**Desired statistical power (\(1 - \beta\)) and false discovery rate (\(\text{FDR}\)).
Socratic Diagnostic Scenarios
Scenario 1: The Log-Transformation Pseudocount Artifact
Question: A researcher log-transforms raw counts (\(\log_2(\text{counts} + 1)\)) and runs a standard Student’s \(t\)-test in Excel across 6 samples. They report thousands of differentially expressed genes, mostly with low read counts (0 to 5 reads). How do you critique this result from first principles?
Diagnostic Answer: Adding a static pseudocount \(+1\) distorts the variance structure of low-count genes. For a gene with \(0\) reads, \(\log_2(0 + 1) = 0\). For a gene with \(5\) reads, \(\log_2(5 + 1) = 2.58\). In low-count regimes, small integer fluctuations (\(0 \to 5\) reads) create artificial multi-fold shifts in log-space. Furthermore, un-normalized raw counts violate \(t\)-test assumptions because library size differences and mean-variance heteroscedasticity are ignored. The reported low-count hits are predominantly artifacts of pseudocount transformation rather than biological regulation.
Scenario 2: Diagnostic Bimodal Dispersion Inspection
Question: You inspect a DESeq2 dispersion plot (Figure 4.2) and observe that raw gene-wise estimates split into two distinct parallel clouds across all mean expression levels, rather than shrinking toward a single trend line. What does this diagnostic signal indicate?
Diagnostic Answer: A bimodal dispersion pattern is a diagnostic clue indicating unmodeled experimental structure—such as an unrecorded batch effect (samples processed on two different dates or sequencing runs) or unexpected subgroup heterogeneity (e.g., male vs. female samples or distinct cell subtype proportions). You must audit your sample metadata, inspect principal components analysis (PCA) plots, and update your design formula (
~ batch + condition) to account for the secondary covariate before finalizing differential expression testing.
The Practical Statistical Audit Checklist (“What Would You Do Tomorrow?”)
When evaluating a differential expression table or designing an omics study, I recommend executing this 7-step statistical audit:
Verify biological replication: Ensure sample groups represent independent biological organisms, not repeated technical measurements of the same sample.
Audit library size normalization: Confirm that raw counts are normalized using size factors (e.g., DESeq2 median ratio or edgeR TMM) to control for sequencing depth.
Inspect the mean-dispersion plot: Verify that dispersion shrinkage was performed and check for anomalous bimodal trend distributions.
Distinguish fold change from statistical evidence: Do not rank genes solely by log fold change or solely by \(p\)-value; evaluate both effect size and uncertainty.
Check dispersion outliers: Ensure high-dispersion genes were not forcibly shrunk down to global trends.
Apply multiple-testing correction: Always adjust raw \(p\)-values using False Discovery Rate (FDR / Benjamini-Hochberg) correction.
Perform power analysis prior to experiment: Use pilot dispersion estimates to calculate required sample sizes for target fold changes.
Current Landscape
Quantitative count modeling continues to advance across single-cell and spatial omics:
Single-Cell Count Modeling (SCTransform v2): Early single-cell RNA-seq algorithms used zero-inflated Poisson or zero-inflated negative binomial models to handle sparse counts. Modern benchmarks (such as
sctransform v2in Seurat) demonstrate that sparse single-cell counts are accurately modeled by regularized Negative Binomial models, using generalized linear models to remove technical sequencing depth variation without zero-inflation parameters.Bayesian Hierarchical Modeling (Stan & brms): Complex multi-factorial omics designs—such as longitudinal clinical trials with patient random effects—increasingly use fully Bayesian hierarchical models. Frameworks built on Stan or PyMC model count distributions while estimating complete posterior uncertainty across biological covariates.
Generalized Additive Models for Trajectories (GAMs & TradeSeq): In single-cell trajectory inference, gene expression changes along continuous pseudotime rather than discrete treatment groups. Tools like
tradeSeqfit Negative Binomial Generalized Additive Models (NB-GAMs), testing differential expression along smooth pseudotime curves.
Summary and Required Reading
Raw counts break un-transformed Gaussian models: Discrete integer counts exhibit mean-variance heteroscedasticity (\(\text{Var}\) scales with \(\mu\)).
Poisson models technical sampling noise: The Poisson constraint (\(\text{Var} = \mu\)) holds for technical replicates of identical libraries.
Biological replicates cause overdispersion: Individual biological variation causes variance to exceed the mean (\(\text{Var} > \mu\)).
Negative Binomial models overdispersion: Characterized by \(\text{Var}(Y) = \mu + \alpha \mu^2\), where \(\alpha\) represents biological dispersion (“room to wander”).
Empirical Bayes MAP shrinkage stabilizes small-\(N\) estimates: DESeq2 fits a mean-dispersion trend \(\alpha(\mu)\) and shrinks gene-wise estimates toward the trend line while preserving high-dispersion outliers.
Effect size is distinct from statistical evidence: A large fold change can be statistically weak if dispersion is high; statistical power depends on effect size, dispersion, depth, and sample size.
Required Reading
- Love, M. I., Huber, W., & Anders, S. (2014): “Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.” Genome Biology, 15(12), 550.
- Robinson, M. D., McCarthy, D. J., & Smyth, G. K. (2010): “edgeR: a Bioconductor package for differential expression analysis of digital gene expression data.” Bioinformatics, 26(1), 139-140.
Johnson’s Rule: Biological replicates are not technical runs, and statistical models are not mathematical conveniences—they are explicit statements about where variation originates. If you model biological sequencing replicates using a Poisson distribution, you are mistaking biological diversity for sampling noise.