10: Quantitative Transcriptomics & Statistical Differential Expression
Chapter 10: Quantitative Transcriptomics & Statistical Differential Expression
RNA sequencing yields discrete integer read counts per gene. However, biological replicates exhibit significant variability exceeding theoretical Poisson expectations. Modeling RNA-seq count data requires specialized generalized linear models (GLMs) governed by the Negative Binomial distribution.
In a standard Poisson distribution, the variance strictly equals the mean (\(\sigma^2 = \mu\)). While technical replicates from the same sequencing library follow Poisson sampling, biological replicates from independent animals or cell cultures exhibit substantial biological variation. The variance systematically exceeds the mean: \(\sigma^2 > \mu\) (overdispersion).
1. The Negative Binomial Model
The Negative Binomial distribution models count variance with a quadratic dispersion parameter \(lpha\):
Var(Kij) = μij + αi × μij2
Where \(\mu_{ij}\) is the expected count for gene \(i\) in sample \(j\), and \(lpha_i\) is the gene-specific dispersion. When sample sizes are small (\(n = 3\) replicates per group), estimating \(lpha_i\) reliably for each gene individually is impossible. Packages like DESeq2 and edgeR use Empirical Bayes shrinkage: borrowing statistical information across all 20,000 genes to shrink gene-wise dispersions toward a common trend line.
2. False Discovery Rate (FDR) & Benjamini-Hochberg
Testing 20,000 genes simultaneously at a standard significance threshold of \(lpha = 0.05\) would yield \(20,000 imes 0.05 = 1,000\) false-positive differentially expressed genes purely by chance. The Benjamini-Hochberg procedure adjusts raw \(p\)-values to control the False Discovery Rate (FDR):
# Benjamini-Hochberg Algorithm:
# 1. Rank all m p-values in ascending order: p(1) <= p(2) <= ... <= p(m)
# 2. Find the largest k such that: p(k) <= (k / m) * alpha
# 3. Reject null hypotheses for all genes ranked 1 through k.
# Adjusted p-value (q-value): q(i) = min [ p(i) * m / i, q(i+1) ]Common Pitfalls & Troubleshooting Checklist
- Java Heap Exhaustion (GATK / Picard): If processes terminate with
java.lang.OutOfMemoryError, set your environment variable explicitly:export _JAVA_OPTIONS="-Xmx16g -Djava.io.tmpdir=/tmp". - Index Timestamp Mismatch: If tools complain that a BAM index is missing or out of date, regenerate with
samtools index -@ 4 sample.bam. - Conda Channel Conflicts: Avoid ABI incompatibilities between bioconda and conda-forge by setting strict priority:
conda config --set channel_priority strict. - Chromosome Naming Contamination: Always verify whether reference FASTA headers use
chr1or1; divergent naming causes silent alignment drops.