We want to show that when \(n\) is large and \(p\) is small such that \(np = \lambda\) remains constant, the \(\text{Binomial}(n, p)\) distribution converges to \(\text{Poisson}(\lambda)\).

Starting with the Binomial PMF:

Substitute \(p = \lambda/n\):

Simplify as \(n \to \infty\):

Therefore:

This approximation is excellent when \(n \geq 100\) and \(p \leq 0.01\), or more generally when \(np < 10\) and \(n\) is large.

We will show that if \(\lambda \sim \text{Gamma}(\alpha, \beta)\) and \(K|\lambda \sim \text{Poisson}(\lambda)\), then marginally \(K \sim \text{NegativeBinomial}\).

Setup:

Marginal distribution by integrating out \(\lambda\):

Recognize the Gamma integral:

Substitute back:

This is the Negative Binomial PMF with:

The variance exceeds the mean by \(\mu^2/\alpha\), capturing over-dispersion.

We show that for \(K \sim \text{NegativeBinomial}\) with mean \(\mu = \alpha/\beta\), the variance is \(\mu + \mu^2/\alpha\).

Using the Gamma-Poisson hierarchy:

Law of Total Expectation:

Law of Total Variance:

For \(\text{Gamma}(\alpha, \beta)\):

Therefore:

When \(\alpha\) is small (high dispersion), the quadratic term \(\mu^2/\alpha\) dominates, creating large over-dispersion. When \(\alpha\) is large, variance approaches \(\mu\) (Poisson limit).

We show that if \(X_1, X_2, \ldots, X_n\) are independent \(\text{Bernoulli}(p)\) random variables, then their sum \(Y = \sum_{i=1}^n X_i\) follows \(\text{Binomial}(n, p)\).

Setup:

Define \(Y = X_1 + X_2 + \cdots + X_n\). We want to find \(P(Y = k)\) for \(k = 0, 1, \ldots, n\).

Counting argument:

\(Y = k\) means exactly \(k\) of the \(n\) trials are successes (value 1) and \(n-k\) are failures (value 0).

Probability calculation:

For any specific sequence with exactly \(k\) successes:

Total probability:

This is exactly the PMF of \(\text{Binomial}(n, p)\). Therefore, the sum of \(n\) independent Bernoulli trials is Binomial.

We prove that for \(X \sim \text{Poisson}(\lambda)\), both \(\mathbb{E}[X] = \lambda\) and \(\text{Var}(X) = \lambda\).

Poisson PMF:

Proof that \(\mathbb{E}[X] = \lambda\):

Proof that \(\text{Var}(X) = \lambda\):

Use \(\text{Var}(X) = \mathbb{E}[X^2] - (\mathbb{E}[X])^2\). First find \(\mathbb{E}[X^2]\):

Therefore:

This unique property-variance equals mean-characterizes the Poisson distribution and is why it's unsuitable for over-dispersed RNA-seq data.

We explain why products of many small random factors lead to a lognormal distribution via the Central Limit Theorem.

Setup: Multiplicative model

Suppose the effective expression rate \(\lambda\) is determined by multiplicative factors:

where each \(F_i > 0\) is a random factor (e.g., library efficiency, cell-type proportion, capture rate).

Taking logarithms:

Central Limit Theorem:

If \(\log(F_i)\) are independent random variables with finite mean \(\mu_i\) and variance \(\sigma_i^2\), then for large \(m\):

Lognormal distribution:

Since \(\log(\lambda)\) is approximately normal, \(\lambda = e^{\log(\lambda)}\) follows a lognormal distribution:

Why Gamma instead?

The lognormal is mathematically valid, but:

• Gamma distributions are conjugate priors for Poisson (analytically tractable)
• Gamma can approximate lognormal well for moderate variation
• The Gamma-Poisson mixture has a closed-form solution (Negative Binomial)
• Lognormal-Poisson mixture has no closed form, requiring numerical integration

Practical takeaway: The multiplicative factors story justifies a heavy-tailed, positive distribution for \(\lambda\). Gamma is chosen for mathematical convenience while capturing the same essential features.

We derive the mean and variance of \(\lambda \sim \text{Gamma}(\alpha, \beta)\).

Gamma PDF:

Mean (First Moment):

Second Moment:

Variance:

Summary:

Note: The variance can be rewritten as \(\text{Var}(\lambda) = \frac{1}{\alpha} \left(\frac{\alpha}{\beta}\right)^2 = \frac{\mu^2}{\alpha}\), showing that smaller \(\alpha\) means higher relative variability.

DESeq2 uses empirical Bayes shrinkage to stabilize dispersion estimates when sample sizes are small.

Problem:

With only 3-6 replicates per condition, the maximum likelihood estimate \(\hat{\theta}_j\) of gene \(j\)'s dispersion is highly variable (high variance estimator).

Observation:

Dispersion tends to depend on mean expression: \(\theta_j = f(\mu_j) + \epsilon_j\) where \(f\) is a smooth function (often power-law or spline).

Empirical Bayes Strategy:

where \(w_j \in [0,1]\) is a data-driven weight that balances:

• Gene-specific evidence: \(\hat{\theta}_j\) (maximum likelihood estimate)
• Prior information: \(f(\mu_j)\) (fitted trend across all genes)

Weight calculation:

The weight \(w_j\) depends on the precision (inverse variance) of \(\hat{\theta}_j\):

where \(\tau^2\) is the variance of true dispersions around the fitted trend (estimated from data).

Effect:

• Low-count genes (high \(\text{Var}(\hat{\theta}_j)\)) → \(w_j \approx 0\) → shrink heavily toward trend
• High-count genes (low \(\text{Var}(\hat{\theta}_j)\)) → \(w_j \approx 1\) → retain gene-specific estimate

Benefit: Reduces false positives from genes with unreliable dispersion estimates while preserving true biological variability.

When testing \(m\) hypotheses simultaneously (e.g., 20,000 genes), we need to control the False Discovery Rate (FDR): the expected proportion of false positives among rejected hypotheses.

Problem with naive thresholding:

If we reject all genes with \(p < 0.05\) and test 20,000 genes where all nulls are true, we expect \(20{,}000 \times 0.05 = 1{,}000\) false positives!

Benjamini-Hochberg (BH) Procedure:

Given \(m\) p-values \(p_1, \ldots, p_m\), to control FDR at level \(\alpha\) (e.g., 0.05):

Step 1: Order p-values: \(p_{(1)} \leq p_{(2)} \leq \cdots \leq p_{(m)}\)

Step 2: Find the largest \(k\) such that:

Step 3: Reject hypotheses \(H_{(1)}, \ldots, H_{(k)}\)

Adjusted p-values:

DESeq2 reports \(\text{padj}_i\), the smallest FDR level at which hypothesis \(i\) would be rejected:

Interpretation:

• \(\text{padj} < 0.05\) means: gene is significant at FDR 5%
• Among all genes with \(\text{padj} < 0.05\), we expect ≤5% to be false positives

FDR Guarantee (under independence or positive dependence):

This is much weaker than family-wise error rate (FWER) control (e.g., Bonferroni) but more powerful for high-throughput screening.

The Mathematical Journey:

Step 1: Express λ as a product

where \(F_i\) represents factors like expression bursts, cell-type proportions, library amplification, etc.

Step 2: Take logarithms (products become sums)

Step 3: Apply Central Limit Theorem

If \(\log(F_i)\) are independent random variables with finite mean \(\mu_i\) and variance \(\sigma_i^2\), then for large \(m\):

Step 4: Exponentiate to get λ

Since \(\log(\lambda)\) is approximately normal, \(\lambda = e^{\log(\lambda)}\) follows a lognormal distribution:

Step 5: Why approximate with Gamma?

• Conjugate prior: Gamma is the conjugate prior for Poisson (analytically tractable)
• Closed-form solution: Gamma-Poisson mixture = Negative Binomial (exact)
• Good approximation: Gamma can approximate lognormal well for moderate variation
• Practical: Lognormal-Poisson mixture has no closed form, requiring numerical integration

Visual comparison:

Both Gamma and Lognormal are positive, right-skewed distributions that can model the same range of shapes. The Gamma wins due to mathematical convenience!

The Key Reason: Mathematical Elegance

Gamma distribution is chosen primarily because it is the conjugate prior for the Poisson distribution.

What does this mean?

The Four Benefits:

1. Closed-form solution (No numerical integration needed)

• Gamma-Poisson → Negative Binomial (exact analytical form)
• Lognormal-Poisson → No closed form (requires numerical approximation)

2. Interpretable parameters

• Shape parameter α: Controls over-dispersion (smaller α = more over-dispersed)
• Rate parameter β: Scales the mean expression level
• Mean: μ = α/β
• Variance: σ² = α/β²

3. Matches biological reality

• Positive support (λ > 0): Can't have negative rates
• Right-skewed: Most genes low λ, few genes high λ
• Flexible shape: Can model different levels of variation

4. Computational efficiency

• Fast likelihood calculations for thousands of genes
• Stable parameter estimation algorithms
• Well-studied statistical properties

Bottom Line:

While lognormal is theoretically justified by multiplicative factors, Gamma is chosen because it leads to the Negative Binomial distribution-giving us a tractable, interpretable, and computationally efficient model for RNA-seq counts!

What is the Wald Test?

The Wald test is a statistical hypothesis test used to determine if a coefficient in our regression model is significantly different from zero.

The Test Statistic:

where:

• \(\hat{\beta}_{j,\text{treatment}}\): Estimated log fold change for gene j
• \(\text{SE}(\hat{\beta}_{j,\text{treatment}})\): Standard error of the estimate

Under the null hypothesis (no differential expression, \(\beta_{j,\text{treatment}} = 0\)):

The test statistic follows a standard normal distribution.

Computing the p-value:

where \(\Phi\) is the standard normal cumulative distribution function. We use a two-tailed test (factor of 2) because we care about both up- and down-regulation.

Interpretation:

• Large |W|: Estimated effect is many standard errors away from zero → Strong evidence against null → Small p-value
• Small |W|: Estimated effect is close to zero relative to uncertainty → Weak evidence → Large p-value

Why Wald Test (instead of Likelihood Ratio Test)?

• Computationally faster (no need to fit reduced model)
• Asymptotically equivalent to LRT under regularity conditions
• Provides both test statistic and standard errors in one step

Important Note:

For genes with very low counts or high dispersion, the asymptotic normal approximation may not be accurate. DESeq2 handles this through filtering and shrinkage.

The Problem: Noisy LFC Estimates for Low-Count Genes

Consider a gene with:

• Control: 2 reads average
• Treatment: 8 reads average
• Estimated log₂(FC) = log₂(8/2) = 2 (4-fold change!)

But is this reliable? With so few counts, this could be noise!

The Solution: Empirical Bayes Shrinkage

DESeq2 shrinks LFC estimates toward zero, with the amount of shrinkage depending on:

1. Standard error: Higher uncertainty → More shrinkage
2. Gene expression level: Lower counts → More shrinkage
3. Dispersion: Higher dispersion → More shrinkage

The Shrinkage Formula:

where:

• \(\hat{\beta}_j\): Original (MLE) log fold change estimate
• \(\text{SE}_j\): Standard error of the estimate
• \(\text{SE}_{\text{prior}}\): Prior standard deviation (estimated from all genes)

Interpretation:

• Low uncertainty (\(\text{SE}_j\) small): Little shrinkage, keep original estimate
• High uncertainty (\(\text{SE}_j\) large): Strong shrinkage toward zero

Example:

Benefits:

• Reduces false discoveries: Prevents calling low-count noise as significant
• Improves ranking: Genes with reliable large effects rise to top
• Better for downstream: More accurate effect sizes for pathway analysis

Important:

Shrinkage affects LFC estimates (effect sizes), not p-values. P-values still come from the original (unshrunken) Wald test.

DESeq2 Function:

Use lfcShrink() after DESeq() to get shrunken LFC for plotting and ranking genes.