The False Discovery Rate

Every hypothesis test can go wrong in two ways: a false positive (calling a null gene significant — a Type I error) and a false negative (missing a truly differential gene — a Type II error). Classical multiple-testing correction, including Bonferroni, targets the family-wise error rate (FWER): the probability that at least one false positive occurs anywhere among all the tests performed. FWER control is appropriate when a single false conclusion would be costly and unacceptable — confirming one specific drug target, for instance. But it scales terribly: as the number of tests m grows into the thousands, the per-test threshold α/m shrinks toward zero, and power to detect real effects collapses. The false discovery rate, introduced by Yoav Benjamini and Yosef Hochberg in 1995, targets a different quantity entirely:

where R is the total number of tests declared significant ("discoveries") and V is the number of those that are actually false positives (true nulls wrongly rejected). FDR is the expected proportion of false positives among your discoveries — not the probability of making any error at all. If you set your procedure to control FDR at 5% and it flags 400 genes, you expect about 20 of them to be false leads and about 380 to be real — you don't know which 20, but you know the rate. This distinction matters enormously in practice. A biologist doing 20,000 tests who accepts a 5% FDR is saying: "I am willing to chase a modest number of false leads through validation experiments, in exchange for catching most of the true signal." A biologist demanding Bonferroni-level FWER control is saying: "I will accept almost no true discoveries in order to guarantee zero false ones." For hypothesis-generating, high-throughput screens — the daily reality of genomics — the FDR philosophy is almost always the right one, because follow-up experiments (qPCR validation, knockout studies, replication cohorts) exist precisely to weed out the false positives that slip through.

The genius of the Benjamini–Hochberg (BH) procedure is that it requires nothing more than the ordinary p-values already computed from each test, plus one simple sorting-and-comparison step.

Step 1 — Compute p-values for all m tests. Run the same statistical test (a t-test, a chi-square test, a likelihood-ratio test) independently on each gene, SNP, or protein, yielding p₁, p₂, …, p_m.

Step 2 — Sort the p-values in ascending order. Relabel them p₍₁₎ ≤ p₍₂₎ ≤ … ≤ p₍m₎, where p₍k₎ is the k-th smallest p-value.

Step 3 — Compute the BH critical value for each rank. For a chosen FDR level q (commonly 0.05 or 0.10), the critical value at rank k is:

This is the key departure from Bonferroni. Bonferroni uses a single flat threshold, q/m, for every test regardless of rank. BH instead uses a linearly increasing staircase of thresholds — the largest p-values must still be very small, but the smallest p-values (rank 1, rank 2 …) are held to almost the same strict standard as Bonferroni, while p-values further down the sorted list are allowed to be progressively larger and still count as significant.

Step 4 — Find the largest rank k where p₍k₎ ≤ (k/m)q. Call this rank k<strong>. All hypotheses with rank ≤ k — that is, all tests with p-values less than or equal to p₍k*₎ — are declared significant.

Step 5 — Reject all null hypotheses at or below that threshold. Note that this is not simply "reject everything with p ≤ (k/m)q applied test-by-test" — it is a single global cutoff, p₍k₎, found by scanning from the largest p-value backward until the staircase condition is first satisfied. Why does this control FDR at level q*? Intuitively: under the null, p-values are uniformly distributed between 0 and 1. If a large fraction of the m tests are true nulls, their p-values scatter roughly evenly across (0,1), so only about q·m of them are expected to fall below any given threshold t by chance, for threshold t = q. The BH staircase exploits this by allowing the threshold to rise proportionally to how many tests you're willing to call significant, so that the expected number of nulls surviving stays proportional to q times the total discoveries — which is exactly the FDR definition. Benjamini and Hochberg proved this holds exactly when the tests are independent (or under certain positive-dependence structures common in real correlated data, as later extended by Benjamini and Yekutieli in 2001).

The figure shows the essential geometry: the green points are p-values that fall below the rising BH staircase — they are declared significant even though several of them sit well above the flat, unforgiving Bonferroni line. Bonferroni, demanding every p-value beat 0.005 regardless of rank, rejects nothing in this toy example; BH, allowing the threshold to climb with rank, recovers eight discoveries while still bounding the expected false-discovery proportion at 10%.

The q-value. Just as a p-value answers "what is the smallest significance level at which this single test would be rejected?", the q-value, introduced by John Storey in 2002, answers "what is the smallest FDR at which this particular test would be included among the discoveries?" Practically, the q-value of the gene ranked k-th is computed by taking p₍k₎·m/k and then enforcing monotonicity (a q-value cannot be smaller than the q-value of a lower-ranked, i.e., smaller, p-value) by taking a running minimum from the largest rank down to the smallest. Reporting q-values lets a reader choose their own acceptable FDR threshold after the fact, rather than being locked into whatever q the original analyst chose — much as reporting a raw p-value lets a reader apply their own alpha.

An important subtlety: FDR control assumes a mixture. Among the m tests, some unknown fraction π₀ are truly null (no real effect) and the rest, 1 − π₀, are truly alternative (real effects). If π₀ is close to 1 — almost nothing in the experiment is real — FDR-controlled discoveries will be few and noisy no matter the method. If π₀ is small — the screen is enriched for real signal, as when comparing tumor vs. normal tissue at genes already suspected of cancer relevance — FDR control is especially powerful, because a large share of small p-values genuinely reflect true effects, letting the staircase climb steeply before hitting a p-value too large to trust.

Benjamini & Hochberg (1995), Journal of the Royal Statistical Society, Series B. The foundational paper, "Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing," did not run a wet-lab biology experiment — it is a statistical proof, but one motivated explicitly by the needs of large-scale applied research (they cite agricultural field trials and clinical multi-endpoint studies as their motivating cases). Benjamini and Hochberg proved, mathematically, that the step-up procedure described above controls E[V/R] at exactly π₀·q when tests are independent, and at most q overall — a guarantee, not a heuristic. Their key theoretical contribution was showing that this control holds regardless of how many of the m null hypotheses are actually true and regardless of the configuration of true and false nulls, which is what makes the procedure broadly applicable rather than dependent on assumptions the experimenter cannot verify in advance.

Benjamini & Yekutieli (2001). Real biological data are not independent — genes on the same pathway are co-regulated, SNPs in linkage disequilibrium are correlated, adjacent brain voxels share signal. Benjamini and Yekutieli extended the 1995 proof to show that the same BH procedure still controls FDR under positive regression dependency (PRDS), a condition satisfied by many realistic correlation structures in genomic data, and that a simple modification (dividing q by the harmonic sum of 1 through m) guarantees control under arbitrary, even negative, dependence — at the cost of some conservatism. This closed the most important practical objection to applying BH to real -omics datasets, where independence is the exception rather than the rule.

Storey (2002) and the q-value. John Storey reframed FDR estimation from a hypothesis-testing procedure into a point-estimation problem: rather than fixing q in advance and asking which tests pass, estimate the proportion of true nulls π₀ directly from the overall p-value distribution (nulls are uniform, so a flat "floor" in the histogram of all p-values estimates π₀), then compute a q-value for every test. This is now the standard output of genomics software (DESeq2, edgeR, limma all report BH-adjusted p-values or q-values alongside raw p-values), letting the same result set support different FDR thresholds for different downstream uses.

Empirical validation in genomics. The clearest demonstration that FDR control behaves as advertised comes from permutation and simulation studies embedded in landmark microarray and RNA-seq papers of the early 2000s (e.q., studies re-analyzing the Golub et al. 1999 leukemia dataset). When the sample labels (tumor vs. normal) are randomly permuted — destroying any true biological signal by construction — applying BH at q = 0.05 to the permuted data yields close to 5% of the "discoveries" as false positives on average across repeated permutations, exactly matching theory. This permutation-based sanity check is now routine practice: it is the empirical proof, repeated in thousands of individual studies, that the mathematical guarantee translates into real behavior on real correlated genomic data.

Example 1 — An RNA-seq differential expression screen. A study compares gene expression in 10,000 genes between drought-stressed and well-watered maize plants. The raw p-values, sorted ascending, include (among 10,000 total) these ten smallest values:

With m = 10,000 and q = 0.05, the threshold at rank k is (k/10000)(0.05) = k × 0.0000050. Checking each rank against its own threshold: Rank 1: 0.0000012 ≤ 0.0000050 ✓ Rank 2: 0.0000031 ≤ 0.0000100 ✓ Rank 3: 0.0000200 ≤ 0.0000150? No — 0.0000200 > 0.0000150. Rank 4 onward: also fail (each p-value is far above its threshold). The BH rule requires finding the largest rank at which the condition holds and rejecting all ranks up to and including it. Ranks 1 and 2 satisfy the condition individually, but rank 3 fails — and since no larger rank in this truncated list satisfies the inequality either (the p-values grow faster than the threshold here), the largest satisfying rank is k\ = 2. <strong>Only the two genes at ranks 1 and 2 are declared significant</strong> at FDR = 5%, even though rank 3's raw p-value (0.00002) would have looked highly significant by ordinary standards. Compare this to Bonferroni: the flat threshold is 0.05/10,000 = 0.000005, which only ranks 1 and 2 pass here too — in this particular (adversarial) data configuration BH and Bonferroni happen to agree, illustrating that BH is never less* conservative than the data structure allows, but in typical genomic datasets with many true small p-values clustered together, BH's staircase pulls far more discoveries through than Bonferroni's flat line (as shown in the first figure).

Example 2 — Choosing an FDR threshold in a GWAS-style hit-picking exercise. A genome-wide association study tests 500,000 SNPs for association with type 2 diabetes risk. The researchers plan two different follow-ups depending on how many candidate SNPs they get: (a) an expensive functional validation in mouse models (can afford to chase ~10 candidates), or (b) a cheap in-silico annotation pipeline (can screen thousands). Suppose the BH procedure at q = 0.05 yields 340 significant SNPs, while at q = 0.01 it yields 38 significant SNPs, and at q = 0.20 it yields 1,150 significant SNPs. (This ordering is guaranteed: lowering q can only shrink or hold constant the discovery set, since the staircase threshold (k/m)q scales down uniformly.) For plan (a), the mouse validation budget of ~10 candidates argues for the strictest threshold available or even a manual scan of the q-value list, picking the top 10 by q-value regardless of which nominal q was used to generate the list — since q-values are monotonic with rank, the top 10 q-values are the same set whether you computed them targeting 0.01 or 0.05. For plan (b), where thousands of candidates can be tolerated and the in-silico pipeline itself has its own filtering power (only truly interesting hits get flagged by annotation), the 5% or even 20% FDR list is defensible: even if roughly 17 of the 340 (5% FDR) hits are false, or roughly 230 of the 1,150 (20% FDR) hits are false, the downstream computational filter can absorb that noise cheaply. This example shows that the "right" FDR threshold is not a universal constant like α = 0.05 — it is a decision-theoretic choice tied to the cost of chasing a false lead versus the cost of missing a true one, which the p-value/q-value machinery makes explicit rather than hiding.

Genome-wide association studies (GWAS). Testing hundreds of thousands to millions of SNPs against a phenotype demands FDR or related control (many GWAS also use a fixed "genome-wide significance" threshold near 5×10⁻⁸, itself derived from Bonferroni-style reasoning over an estimated number of independent tests, but FDR/q-values are widely reported alongside it for ranking candidate loci for replication).

RNA-seq differential expression. Standard tools — DESeq2, edgeR, limma-voom — report BH-adjusted p-values by default for every gene in a transcriptome-wide comparison, letting researchers set an FDR cutoff (typically 0.05 or 0.10) to generate a manageable, mostly-valid gene list for pathway analysis.

Proteomics and metabolomics. Mass-spectrometry-based peptide identification screens tens of thousands of candidate peptide-spectrum matches; FDR control (often via a "target-decoy" strategy, searching a reversed or scrambled database to empirically estimate the false-match rate) is the industry standard for reporting confident protein identifications.

Epigenomics and ChIP-seq. Peak-calling software (e.g., MACS2) reports q-values for each candidate binding site or histone-modification peak across the genome, since millions of genomic bins are implicitly tested for enrichment.

Neuroimaging. fMRI studies testing tens of thousands of brain voxels for activation differences between conditions routinely use FDR-based thresholds (alongside cluster-based and permutation methods) to balance sensitivity against the "voodoo correlations" and false-positive scandals (e.g., the 2009 "dead salmon" fMRI study) that motivated stricter multiple-testing discipline in the field.

Clinical trials with multiple endpoints. Trials measuring several secondary outcomes (biomarkers, sub-scale scores) increasingly report FDR-adjusted p-values for exploratory endpoints, reserving strict FWER control for the pre-specified primary endpoint.

Conservation and ecological screens. Comparative genomics scans for signatures of selection across thousands of loci in endangered or invasive species use FDR control to prioritize a tractable number of candidate loci for functional follow-up, rather than either drowning in false Bonferroni-purged emptiness or chasing every marginally elevated statistic. Across all these domains, FDR control is what makes exploratory, hypothesis-generating "-omics" science statistically responsible without being paralyzingly conservative — a bridge between the rigor demanded by the Multiple-Testing Problem and the practical need to actually find something, setting the stage for the broader statistical machinery used in Statistics of High-Dimensional Data and in Genome-Wide Association Statistics.

Ask the tutor: Why does the Benjamini–Hochberg procedure remain valid even when I don't know how many of my tested genes are truly null?