What the Benjamini-Hochberg Procedure Does, and How to Use It on Your Own Omics Data

Benjamini-Hochberg (BH) is a procedure that takes a list of p-values from many hypothesis tests and decides which to call significant, such that the expected fraction of false positives among your calls stays below a level you choose. That level is the false discovery rate, or FDR. FDR is the error rate, and BH is one method for controlling it. BH and FDR are not the same thing. BH is the method almost every RNA-seq and proteomics pipeline uses by default, and it is the right default for most personal omics work.

The procedure, step by step

You have $m$ p-values and a target rate $q$, commonly 0.05 or 0.10. Sort the p-values in ascending order so that $p_{(1)} \le p_{(2)} \le \dots \le p_{(m)}$. Find the largest index $k$ for which $p_{(k)} \le \frac{k}{m} q$. Reject the null hypothesis for every test from rank 1 through rank $k$, including any whose own p-value exceeds its threshold.

That last clause is what makes BH a step-up procedure rather than a simple per-test cutoff. The threshold slides: the smallest p-value must beat $q/m$, which is the Bonferroni bar, while the largest must only beat $q$. Software usually reports BH-adjusted p-values instead of a yes-or-no call. The adjusted value at rank $i$ is $\min_{j \ge i} \left( \frac{m}{j} p_{(j)} \right)$, a running minimum taken from the largest p-value downward. This enforces monotonicity so adjusted values never decrease as raw p-values increase. Thresholding the adjusted values at $q$ gives exactly the same rejection set as the procedure above.

BH controls FDR at $q$ under independence of the test statistics, and also under a dependence condition called positive regression dependency. Gene expression is heavily correlated, so strict independence never holds. The practical situation is better than that sounds, because positive correlation among co-regulated features usually falls inside the class BH tolerates. Arbitrary dependence requires the Benjamini-Yekutieli variant, which divides $q$ by $\sum_{i=1}^{m} 1/i$ and costs roughly a factor of 7 in threshold at $m = 1000$. We rarely use it for expression data, and the reasons for and against it are laid out well in recent reviews of FDR methods beyond BH.¹

Bonferroni, and why FDR replaced it for omics

Bonferroni answers a different question. It controls the family-wise error rate: the probability of making even one false positive across the whole experiment. You get it by comparing each p-value to $\alpha/m$, or equivalently by multiplying each p-value by $m$ and capping at 1. With 20 tests and $\alpha = 0.05$ each test must clear 0.0025. That is reasonable. With 20,000 genes each test must clear $2.5 \times 10^{-6}$, which in a study of a handful of samples means you will find nothing.

The tradeoff is explicit. Bonferroni protects a claim that your entire result set is clean. BH protects a claim about a proportion: at $q = 0.10$, roughly 10% of the features on your list are expected to be false. For screening experiments that generate hypotheses rather than settle them, the second guarantee is the useful one. Simulation work on realistic omics models shows how sharply sensitivity differs among correction methods at fixed effect size and sample count.² Follow-up work extends that comparison across three models of quantitative omics experiments and quantifies the power lost to family-wise control.³

Running it

In R, the one-liner is p.adjust(pvals, method = "BH"). Note that p.adjust uses the length of the input vector as $m$, so any filtering you intend must happen before the call. method = "BY" gives Benjamini-Yekutieli and method = "bonferroni" gives Bonferroni.

In Python, use statsmodels.stats.multitest.multipletests(pvals, alpha=0.10, method="fdr_bh"), which returns a boolean rejection array, the adjusted p-values, and two corrected alpha levels. Pass method="fdr_by" for the arbitrary-dependence version. If you are working from DESeq2 output, the padj column is already BH-adjusted and the pvalue column is not. In limma, topTable(fit, adjust.method = "BH", number = Inf) gives the same in the adj.P.Val column.

Two adjacent tools are worth knowing. Storey’s q-value estimates $\pi_0$, the proportion of true nulls, and divides the BH threshold by it. When a large share of features are genuinely changed, this recovers real power, and qvalue::qvalue() in R computes it. Independent hypothesis weighting (IHW) upweights tests by an independent covariate such as mean expression, and IHW::ihw() typically yields 10–30% more discoveries on RNA-seq at the same nominal FDR. When your features fall into natural groups, group-aware variants exist, including a two-dimensional procedure for two-way designs.⁴

The three failure modes that matter

The first is miscalibrated p-values. BH assumes that null p-values are uniform on the interval from 0 to 1. Plot a histogram of your raw p-values with 50 bins before you adjust anything. You want a flat body with a spike in the leftmost bin. A hump in the middle means the null model is wrong, often from unmodeled batch structure or overdispersion. A slope rising toward 1 means the test is conservative. In either case, no correction can rescue the output, and the fix is the model.

The second is choosing $m$ badly. Every hopeless test you include dilutes the threshold for the tests that could have worked. Filtering out features with near-zero counts or proteins quantified in fewer than half the samples is legitimate and increases power. The filter must not depend on the test statistic, which is why DESeq2’s independent filtering keys off mean normalized count rather than the p-value. Filtering on fold change first, then running BH, invalidates the guarantee.

The third is running BH more than once on overlapping families. Suppose you test 20,000 transcripts and then separately test 7,000 proteins, adjusting each set at $q = 0.05$. You have controlled FDR within each assay and not across the combined claim. Multi-omics studies handle this by declaring the families in advance and reporting per-omic adjusted values. Nested case-control plasma biomarker work in pancreatic cancer illustrates the reporting convention.⁵ Proteomics and metabolomics carry their own wrinkles from missingness and correlated features, which the mass-spectrometry multiple-comparison literature treats directly.⁶

What BH does not tell you

An FDR value is a property of a list. It does not belong to an individual feature. A value of $q = 0.04$ for a gene means that if you accept every gene down to that rank, about 4% of them are expected to be false. The gene is not 96% likely to be real. Two other things follow. Effect size is not addressed at all: a transcript can pass at $q = 10^{-8}$ with a log2 fold change of 0.05, which is measurable and biologically inert. And BH says nothing about whether your comparison answers a question you care about, which is where general principles for omics analysis are more useful than any single correction.⁷

For your own longitudinal data, the framing is that $n = 1$ time series give you within-person contrasts with real power at repeated sampling, and exploratory status otherwise. Exploratory multi-omics profiling in small cohorts is routinely reported that way.⁸ Anything that touches clinical decisions belongs with a physician or a board-certified genetic counselor, and germline variant interpretation uses ACMG criteria rather than any FDR procedure.

Questions people also ask

What is the Benjamini-Hochberg method? A step-up procedure that sorts $m$ p-values, compares the $k$-th smallest to $\frac{k}{m} q$, and rejects everything up to the largest $k$ that passes. It controls the expected proportion of false positives among rejected hypotheses at $q$.

Is Benjamini-Hochberg the same as FDR? No. The false discovery rate is the quantity being controlled. Benjamini-Hochberg is the most common procedure for controlling it, and alternatives include Storey’s q-value, Benjamini-Yekutieli, IHW, and knockoff-based methods.

How do I cite the Benjamini-Hochberg method? Cite Y. Benjamini and Y. Hochberg (1995), “Controlling the false discovery rate: a practical and powerful approach to multiple testing”. Journal of the Royal Statistical Society Series B, 57(1), 289-300.

Can you explain Bonferroni simply? Divide your significance level by the number of tests, or equivalently multiply each p-value by that number. It bounds the chance of any false positive anywhere in the experiment, and at omics scale it is usually too strict to find anything.

Which should I use for RNA-seq? BH at $q = 0.05$ or $q = 0.10$, applied after count-based independent filtering, with a minimum fold-change threshold applied separately for interpretation. Reserve Bonferroni for small confirmatory sets of pre-specified tests.

Footnotes

1. Salil Koner, Navonil Sarkar, Nilanjana Laha. False discovery rate control: Moving beyond the Benjamini–Hochberg method. 2024. https://doi.org/10.21203/rs.3.rs-3861673/v1 ↩

2. Steven R. Shuken, Margaret W. McNerney. Visualizing the Costs and Benefits of Correcting P-Values for Multiple Hypothesis Testing in Omics Data. 2021. https://doi.org/10.1101/2021.09.09.459558 ↩

3. Steven R. Shuken, M. Windy McNerney. Costs and Benefits of Popular P-Value Correction Methods in Three Models of Quantitative Omic Experiments. Analytical Chemistry, 2023. https://doi.org/10.1021/acs.analchem.2c03719 ↩

4. Lu Yang, Pei Wang, Jun Chen. 2dGBH: Two-dimensional group Benjamini–Hochberg procedure for false discovery rate control in two-way multiple testing of genomic data. Bioinformatics, 2024. https://doi.org/10.1093/bioinformatics/btae035 ↩

5. Emmy Borgmästars, Benjamin Ulfenborg, Mattias Johansson, et al. Multi-omics profiling to identify early plasma biomarkers in pre-diagnostic pancreatic ductal adenocarcinoma: a nested case-control study. Translational Oncology, 2024. https://doi.org/10.1016/j.tranon.2024.102059 ↩

6. Pietro Franceschi, Marco Giordan, Ron Wehrens. Multiple comparisons in mass-spectrometry-based -omics technologies. TrAC Trends in Analytical Chemistry, 2013. https://doi.org/10.1016/j.trac.2013.04.011 ↩

7. Daniela Dunkler, Fátima Sánchez-Cabo, Georg Heinze. Statistical Analysis Principles for Omics Data. Methods in Molecular Biology, 2011. https://doi.org/10.1007/978-1-61779-027-0_5 ↩

8. Na Chen, Yubing Hu, Tianxia Xiao, et al. Integrative serum proteomic and metabolomic profiling in ovarian endometrioma: an exploratory multi-omics study. Frontiers in Medicine, 2026. https://doi.org/10.3389/fmed.2026.1846827 ↩