Why Gene-Wise Tests Need a Multiple-Testing Plan
A pseudobulk differential-expression analysis may test thousands of genes for the same cell type and condition contrast. Even if every null hypothesis is true and each gene-wise test is calibrated, some small p-values will occur by chance. The more tests conducted, the more likely it is that an unadjusted threshold such as 0.05 will label at least one gene a discovery.
The multiple-testing procedure addresses this problem across a defined family of hypotheses. It does not replace the gene-wise model or repair a poorly calibrated test. The p-values supplied to the procedure must come from appropriate pseudobulk tests, with the experimental units, design, and contrast specified correctly. Here, the focus is how to adjust those p-values and interpret the results.
FDR control at level \(\alpha\) means that, under the procedure’s assumptions, the expected false discovery proportion is at most \(\alpha\). It does not mean that each selected gene has probability \(\alpha\) of being false, nor does it identify which particular discoveries are false. It is also distinct from family-wise error control, which limits the probability of making one or more false rejections.
Benjamini–Hochberg Adjustment
The Benjamini–Hochberg (BH) procedure is a standard way to control FDR across a family of gene-wise tests. Let \(m\) be the number of tests in that family and let their p-values, sorted from smallest to largest, be \(p_{(1)}\leq p_{(2)}\leq\cdots\leq p_{(m)}\). For a chosen FDR level \(\alpha\), compare each ordered p-value with its rank-specific threshold.
Find the largest rank \(k\) for which the inequality holds. Reject the hypotheses with ranks \(1\) through \(k\); if no rank passes, reject none. The rule is a step-up procedure: a passing rank permits rejection of all hypotheses ranked below it, even if a lower-ranked comparison would not pass when considered in isolation.
The same procedure can be reported through BH-adjusted p-values, often called q-values in software output. For the ordered p-values, calculate \(m p_{(i)}/i\), then enforce nondecreasing adjusted values from the smallest to the largest p-value by taking a cumulative minimum from the largest rank backward. Cap values at 1. A hypothesis is rejected at level \(\alpha\) when its adjusted p-value is at most \(\alpha\).
The cumulative minimum matters: the raw rank-scaled quantities \(m p_{(i)}/i\) need not be ordered, whereas BH-adjusted p-values are nondecreasing with the ordered raw p-values. These adjusted values summarize the BH rejection decisions across possible FDR thresholds; they are not posterior probabilities that individual null hypotheses are true.
Worked Example: Finding the BH Rejection Set
Suppose ten genes in one hypothetical pseudobulk contrast have ordered p-values \(0.004, 0.012, 0.028, 0.041, 0.062, 0.110, 0.200, 0.310, 0.420,\) and \(0.600\). Assume these are valid gene-wise p-values for the stated family. Apply BH at \(\alpha=0.10\).
Here \(m=10\), so the rank-specific thresholds \(i\alpha/m\) are \(0.010, 0.020, 0.030, 0.040, 0.050, 0.060, 0.070, 0.080, 0.090,\) and \(0.100\). At ranks one to four, the comparisons are \(0.004\leq0.010\), \(0.012\leq0.020\), \(0.028\leq0.030\), and \(0.041>0.040\). The first three ranks pass, including rank three; rank four fails. Every later p-value also exceeds its threshold, so the largest passing rank is \(k=3\). BH rejects the first three hypotheses.
For example, the raw rank-scaled value at rank three is \(10(0.028)/3=0.0933\). At rank four it is \(10(0.041)/4=0.1025\). Taking the cumulative minimum from the largest rank backward gives adjusted values for ranks one through four of \(0.040, 0.060, 0.0933,\) and \(0.1025\). The first three are at most 0.10; the fourth is not. The conclusion agrees with the step-up rule.
Why the BH Rule Works
A useful justification shows how ordering and the rank threshold limit the expected false discovery proportion. Under independent tests with valid null p-values, imagine there are \(m_0\) true null hypotheses among \(m\) tests. BH’s expected false discovery proportion is bounded by \(\alpha m_0/m\), which is no larger than \(\alpha\). The bound accounts for the fact that a true-null p-value must be small enough to fall among the rejected ranks, rather than simply comparing every p-value with the same unadjusted cutoff.
The key intuition is that a rejection at rank \(i\) requires a p-value no larger than \(i\alpha/m\). Under a true null, validity gives \(P(p\leq t)\leq t\); independence allows the expected contribution of each true-null test to be bounded while accounting for the rank that test occupies. Summing those contributions yields the stated bound. Standard positive-dependence conditions also support BH control, but arbitrary dependence does not follow from this independent-test justification. For arbitrary dependence, a conservative option is the Benjamini–Yekutieli threshold \(i\alpha/(m H_m)\), where \(H_m=\sum_{j=1}^m 1/j\); it makes the rank thresholds smaller.
Worked Calculations and Reporting
Worked Example: Calculating Adjusted P-Values
Consider six hypothetical gene-wise p-values, already sorted: \(0.003, 0.018, 0.045, 0.120, 0.310,\) and \(0.600\). The rank-scaled quantities \(m p_{(i)}/i\), with \(m=6\), are:
\(6(0.003)/1=0.018\); \(6(0.018)/2=0.054\); \(6(0.045)/3=0.090\); \(6(0.120)/4=0.180\); \(6(0.310)/5=0.372\); and \(6(0.600)/6=0.600\). These are already nondecreasing, so they are also the BH-adjusted p-values.
At \(\alpha=0.05\), only the first adjusted p-value is at most 0.05. Thus only the gene with raw p-value 0.003 is rejected for this six-test family. The final rank-scaled value is \(0.600\), since \(6(0.600)/6=0.600\), not a value larger than the original p-value.
Worked Example: Reporting an Effect With Its Uncertainty
Suppose a hypothetical gene has an estimated contrast of \(0.80\) log2-fold-change units and standard error \(0.25\) log2-fold-change units from a pseudobulk model. Assume the gene-wise Wald test is appropriate and its large-sample normal approximation is adequate. Its Wald statistic is \(z=0.80/0.25=3.20\), giving a two-sided p-value of approximately \(0.0014\). Suppose the gene is part of a five-gene testing family whose ordered p-values are \(0.0014, 0.012, 0.041, 0.200,\) and \(0.600\).
The BH rank-scaled values are \(5(0.0014)/1=0.0070\), \(5(0.012)/2=0.030\), \(5(0.041)/3=0.0683\), \(5(0.200)/4=0.250\), and \(5(0.600)/5=0.600\). They are nondecreasing, so the gene’s BH-adjusted p-value is 0.0070. It is selected at FDR level 0.05.
A conventional 95% Wald confidence interval for the gene’s estimate is \(0.80\pm1.96(0.25)=[0.31,1.29]\) log2-fold-change units. On the fold-change scale, this is approximately \([2^{0.31},2^{1.29}]=[1.24,2.45]\). The interval describes uncertainty in the estimated contrast, while the adjusted p-value concerns multiplicity across the five tests. This is an ordinary gene-wise interval, not a simultaneous interval adjusted to cover all five effects together. A clear report gives the effect estimate, its interval, the raw p-value, the BH-adjusted p-value, and the family and FDR level.
Choosing the Family, Checking Assumptions, and Comparing Methods
The family must match the inferential claim. If the claim concerns genes for one contrast in one cell type, a natural family is all genes tested for that contrast and cell type. If conclusions will be made across multiple cell types or contrasts, adjusting separately within each subset does not control FDR across the combined set. The analyst should decide and report the family in light of the intended claims, rather than choosing a boundary after seeing which genes are significant.
Filtering can reduce the family size and improve power, but it should not select genes because their observed condition differences or p-values are favorable. Filtering based on information independent of the test statistic under the null—such as a prespecified expression threshold—is generally safer than filtering on the same evidence used for testing. The tested-gene count \(m\) must be the number of hypotheses actually included in the declared family.
| Procedure | Target error rate | Typical consequence |
|---|---|---|
| Benjamini–Hochberg | Expected false discovery proportion | Usually permits more discoveries than family-wise control; relies on valid p-values and suitable dependence conditions |
| Bonferroni | Probability of at least one false rejection in the family | Reject when \(p\leq\alpha/m\); valid under arbitrary dependence for valid p-values, but often conservative |
| Benjamini–Yekutieli | Expected false discovery proportion under arbitrary dependence | Uses smaller thresholds \(i\alpha/(mH_m)\); can be substantially more conservative |
Bonferroni and BH answer different questions. For Bonferroni, dividing \(\alpha\) by \(m\) ensures that the probability of at least one false rejection is at most \(\alpha\), by the union bound. This can be appropriate when any false positive has a high cost. BH instead targets the average proportion of false discoveries among selections, often a more relevant goal for exploratory gene screening. Neither procedure makes invalid gene-wise p-values valid.
Diagnostics begin with the gene-wise tests: check model fit, dispersion behavior, design specification, and whether the test’s reference distribution is plausible. Then check the multiplicity bookkeeping: confirm the family, the number of tested genes, the filtering rule, and the adjustment method. Dependence among gene tests is a limitation to consider; do not assume that a large number of genes makes the tests independent. If the family or adjustment was chosen after inspecting results, the reported error-rate interpretation may not hold as intended.
A small adjusted p-value is not a measure of biological importance. A modest effect can be highly significant with precise estimation, while a larger estimated effect can remain uncertain when data are sparse. Conversely, an adjusted p-value above the threshold does not establish that the gene has no effect. Report estimated effects and their uncertainty alongside adjusted p-values, and distinguish statistical selection from scientific interpretation.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- For eight hypotheses tested at \(\alpha=0.05\), what is the BH threshold at rank four?
- With ten ordered p-values, suppose ranks one through three pass their BH thresholds, rank four fails, and all later ranks fail. Which hypotheses are rejected?
- In a six-test family, the largest ordered p-value is 0.600. What is its raw rank-scaled BH value at rank six?
- What is the difference between controlling FDR and controlling family-wise error?
- Why should a gene-wise confidence interval usually be identified as unadjusted when reported beside BH-adjusted p-values?
Solutions. (1) The threshold is \(4(0.05)/8=0.025\). (2) BH is step-up: reject the hypotheses at ranks one through three, since rank three is the largest passing rank. (3) The value is \(6(0.600)/6=0.600\). (4) FDR bounds the expected proportion of false rejections among all rejections; family-wise error bounds the probability of at least one false rejection. (5) BH adjusts the p-values for a family of tests, but does not by itself make the separate gene-wise confidence intervals simultaneous intervals with family-wide coverage.