Differential gene expression (DGE) analysis identifies genes whose expression levels differ significantly between experimental conditions (e.g., treated vs. control, diseased vs. healthy). Tools like DESeq2 and edgeR model RNA-seq count data using the negative binomial distribution (which accounts for both sampling noise and biological variability), estimate per-gene dispersion, and perform statistical tests for each gene. Because thousands of genes are tested simultaneously, multiple testing correction (Benjamini-Hochberg FDR) is essential to control the false discovery rate. Results are typically reported as log2 fold changes with adjusted p-values and visualized using volcano plots and MA plots.
Run DESeq2 on a small RNA-seq dataset with 3 replicates per condition. Examine the results table: sort by adjusted p-value, filter by fold change, and generate a volcano plot. Then repeat the analysis removing one replicate per condition and observe how statistical power decreases — this demonstrates why biological replication matters more than sequencing depth.
The RNA-seq pipeline produces a matrix of read counts per gene per sample. The next question — which genes are expressed differently between conditions? — is fundamentally a statistical problem. Differential gene expression analysis applies statistical models to this count data to identify genes with expression changes larger than expected from random variation.
The statistical framework starts with the right distribution for count data. RNA-seq counts are not normally distributed; they are discrete, non-negative, and often skewed. The Poisson distribution is a natural starting point (it models count data), but it assumes the variance equals the mean. In practice, biological replicates show more variability than Poisson predicts — called overdispersion. The negative binomial distribution adds a dispersion parameter to capture this biological variability. DESeq2 and edgeR both use the negative binomial model but differ in exactly how they estimate dispersion and normalize data.
A critical technical challenge is dispersion estimation with few replicates. With only 3 replicates per condition (common in RNA-seq), estimating the variance of each gene independently would be very noisy. Both DESeq2 and edgeR address this by "borrowing strength" across genes: they assume genes with similar expression levels have similar dispersion, fitting a trend of dispersion versus mean expression and shrinking individual gene estimates toward this trend. This approach — empirical Bayes shrinkage — stabilizes variance estimates and improves statistical power, but it assumes the dispersion-mean relationship is smooth, which generally holds in practice.
Multiple testing correction is non-negotiable. Testing 20,000 genes means that even with well-calibrated p-values, a 5% significance threshold produces ~1,000 false positives. The Benjamini-Hochberg procedure converts p-values to adjusted p-values (q-values) that control the false discovery rate — the expected proportion of false positives among all declared significant results. An FDR threshold of 0.05 means you accept that approximately 5% of your significant genes may be false discoveries, which is a reasonable tradeoff in exploratory genomics where downstream validation (qPCR, functional assays) will filter the list further.
Results are typically visualized with volcano plots (log2 fold change on x-axis, -log10 adjusted p-value on y-axis), which simultaneously show effect size and statistical significance, making it easy to identify genes that are both biologically meaningful (large fold change) and statistically reliable (low adjusted p-value). The output gene list feeds into pathway analysis, gene ontology enrichment, and network analysis to interpret the biological significance of expression changes.