Why Gene-Level Significance Needs Context
An RNA-seq differential-expression analysis can test thousands of genes for a condition effect. Even if every null hypothesis is true, some small p-values can occur by chance. A gene-level p-value describes evidence for that gene under its test, but does not by itself account for the many other genes tested. Adjusted p-values address this multiplicity; they do not tell us whether a change is large, precisely estimated, or biologically consequential.
The relevant interpretation therefore brings together three pieces of information: the adjusted p-value, the estimated effect size, and uncertainty in that estimate. The adjusted p-value concerns evidence across a specified family of tests. The effect size describes the estimated direction and magnitude of change. An interval describes the range of effects reasonably compatible with the data under the model and interval procedure. These answer different questions and should not be substituted for one another.
What the Benjamini–Hochberg Adjustment Does
Let \(m\) be the number of hypotheses in the family being adjusted, and let \(p_{(1)} \le \cdots \le p_{(m)}\) be their ordered, valid p-values. The Benjamini–Hochberg (BH) procedure controls the false discovery rate: the expected proportion of false rejections among all rejections, with that proportion defined as zero when there are no rejections. Under independent tests, or certain forms of positive dependence, the procedure controls this rate at the chosen level \(\alpha\), provided the null p-values are valid (at least as large as a Uniform(0,1) variable in distribution).
At level \(\alpha\), find the largest rank \(k\) for which \(p_{(k)} \le k\alpha/m\). Reject the hypotheses with ranks 1 through \(k\). The corresponding BH adjusted p-value for the hypothesis at rank \(i\) is
The minimum over ranks at or above \(i\) enforces the ordering needed for the step-up procedure. Capping at 1 keeps the result within the p-value range. Tied p-values can be handled by assigning their shared value the appropriate rank-based results; software generally implements a consistent tie convention.
Why the Adjustment Formula Has a Minimum
The formula follows directly from the step-up rejection rule. For a fixed gene at rank \(i\) to be included among the rejections, a cutoff rank \(k\) at least as large as \(i\) must satisfy \(p_{(k)} \le k\alpha/m\). Rearranging gives \(\alpha \ge m p_{(k)}/k\). The smallest level \(\alpha\) at which this can happen is the smallest of these quantities over \(k \ge i\). If that smallest level exceeds 1, it is reported as 1. This yields the BH adjusted p-value above.
This interpretation explains why the adjusted p-value is not simply \(m\) times the raw p-value: the rank and the subsequent p-values matter. It also explains why adjusted values are made nondecreasing with rank, even if the uncorrected quantities \(m p_{(i)}/i\) fluctuate.
Worked Example: Calculating BH Adjusted P-Values
Suppose six genes have the following hypothetical p-values, already sorted from smallest to largest: 0.001, 0.012, 0.018, 0.031, 0.20, and 0.70. The family is these six tests, so \(m=6\).
For rank 1, the unmonotonized value is \(6(0.001)/1=0.006\). For rank 2 it is \(6(0.012)/2=0.036\); for rank 3, \(6(0.018)/3=0.036\); for rank 4, \(6(0.031)/4=0.0465\); for rank 5, \(6(0.20)/5=0.24\); and for rank 6, \(6(0.70)/6=0.70\). Taking the minimum over this and all later ranks gives adjusted values 0.006, 0.036, 0.036, 0.0465, 0.24, and 0.70.
At a target FDR level of \(\alpha=0.05\), the first four adjusted p-values are at most 0.05, while the last two are not. Thus the BH procedure rejects the first four hypotheses at that target level. The fourth gene’s adjusted value is 0.0465, not its raw p-value of 0.031: the adjusted value reflects the family of six tests.
Read Adjusted P-Values With Effect Sizes and Intervals
For a two-condition RNA-seq comparison, an estimated log2 fold change \(\hat{\beta}\) describes the estimated change on a base-2 logarithmic scale. A value of 0 corresponds to no estimated change; a value of 1 corresponds to a twofold increase, and a value of -1 to a twofold decrease. Transforming by \(2^{\hat{\beta}}\) gives the estimated fold change. For example, a log2 fold change of 0.1 corresponds to \(2^{0.1}\), or about a 1.072-fold increase.
An interval for the log2 fold change communicates estimation uncertainty on that scale. Transforming both endpoints by \(2^x\) gives an interval for the fold change. A conventional 95% interval has a repeated-sampling interpretation under its assumptions: in repeated applications of the procedure, 95% of such intervals would cover the fixed parameter. It does not mean there is a 95% probability that this particular interval contains the parameter. Unless constructed with a multiplicity adjustment, a 95% interval is ordinarily a per-gene interval, not a simultaneous 95% interval for all genes.
Worked Example: A Small but Precisely Estimated Change
Suppose a hypothetical gene has estimated log2 fold change \(\hat{\beta}=0.015\), standard error \(SE=0.003\), and a BH adjusted p-value of 0.00004 in a large gene family. The reported adjusted p-value is taken as a hypothetical result from the complete family of tests; it cannot be reconstructed from this gene’s p-value alone.
For illustration, a Wald statistic is \(z=\hat{\beta}/SE=0.015/0.003=5\). A two-sided normal-reference p-value for \(z=5\) is about \(5.73\times10^{-7}\), consistent with strong evidence against a zero coefficient. The adjusted p-value is also small, so the gene passes a 0.05 BH threshold.
A nominal 95% Wald interval on the log2 scale is \(0.015\pm1.96(0.003)\), or \(0.00912\) to \(0.02088\). On the fold-change scale this is \(2^{0.00912}\) to \(2^{0.02088}\), approximately 1.006 to 1.015. The estimated increase is \(2^{0.015}-1\), about 1.05%; the interval corresponds to increases of roughly 0.6% to 1.5%. If a 20% change had been specified in advance as the minimum biologically important increase, its fold change of 1.2 is well above this interval. The small adjusted p-value indicates statistical evidence for a nonzero change, not evidence for a change of that biological size.
Worked Example: Interpreting a Larger Estimated Change
Consider a second hypothetical gene with \(\hat{\beta}=0.12\) log2 fold change and \(SE=0.035\). Suppose its BH adjusted p-value is 0.012. For the coefficient test, the Wald statistic is \(z=0.12/0.035=3.429\), giving a two-sided normal-reference p-value of approximately 0.00061. The adjusted value is an output of correcting the full family of gene tests; here it indicates that the gene passes a 0.05 BH threshold.
A nominal 95% interval is \(0.12\pm1.96(0.035)\). Since \(1.96(0.035)=0.0686\), the interval is 0.0514 to 0.1886 on the log2 scale. Transforming the endpoints gives \(2^{0.0514}\approx1.036\) and \(2^{0.1886}\approx1.140\). Thus, under the model and interval procedure, the interval corresponds to an increase of approximately 3.6% to 14.0%, with a point estimate of \(2^{0.12}\approx1.087\), or about an 8.7% increase.
The adjusted p-value supports evidence against a zero effect across the tested family, while the estimated magnitude is modest. Whether that size matters depends on the scientific question and a biologically meaningful threshold, not on the adjusted p-value alone. The interval is also nominal for this gene unless an interval procedure specifically accounts for multiplicity.
Interpreting Thresholds and the Tested Family
An adjusted p-value of 0.04 means that this gene would be rejected by the BH procedure at FDR levels of 0.04 or higher, and not at lower levels, subject to the procedure’s assumptions and the family used. It does not mean there is a 4% probability that the gene is a false discovery. Nor does an FDR target of 0.05 mean that exactly 5% of the rejected genes in this particular analysis are false. The FDR is an expectation over repeated applications of the procedure, not a guarantee about the realized fraction in every data set.
The family must be specified: for example, all tested genes for a particular contrast, or a broader collection of tests if that is the scientific error-control goal. Adjusted values from different contrasts or analyses are not automatically comparable if their families or filtering procedures differ. Independent filtering can change which genes enter the adjusted family; the interpretation should reflect the analysis that produced the reported values.
Worked Example: A Result Just Above the Chosen Threshold
Suppose four hypothetical genes in one family have sorted p-values 0.004, 0.018, 0.041, and 0.20. Here \(m=4\). The BH quantities are \(4(0.004)/1=0.016\), \(4(0.018)/2=0.036\), \(4(0.041)/3\approx0.0547\), and \(4(0.20)/4=0.20\). They are already nondecreasing, so these are the adjusted values.
At \(\alpha=0.05\), the first two genes pass, while the third does not because \(0.0547>0.05\). At \(\alpha=0.10\), the third passes. Its adjusted p-value near 0.055 is therefore evidence that the BH rejection threshold would need to be at least that large for it to be selected. It is not evidence that the null hypothesis has a 5.5% probability of being true. Whether the third gene is scientifically interesting still requires its effect estimate, uncertainty, and biological context.
Diagnostics, Alternatives, and Limitations
Adjustment is a final inference step that relies on the gene-level p-values being appropriate for the tested hypotheses. It cannot repair an unsuitable count model, an incorrectly specified design, unaddressed confounding, or influential observations that invalidate the test calibration. The residual and count diagnostics used earlier in the analysis remain relevant. A large collection of apparently tiny p-values, or a p-value distribution with an unexpected shape, can prompt review of model assumptions, sample structure, and filtering rather than merely a more aggressive adjustment.
Check what the software adjusted, how many genes were in the family, whether filtering was applied, and which procedure was used. A report containing only the adjusted p-value is incomplete for interpretation: include the contrast, estimated log2 fold change (or fold change), and an uncertainty measure. Check whether reported intervals are per-gene or multiplicity-adjusted, and do not pair an interval for one estimator with a different effect estimate without explaining the distinction. Shrinkage can stabilize effect estimates, but it does not convert an adjusted p-value into a measure of effect size.
BH is commonly used because it controls FDR rather than the probability of making even one false rejection. Bonferroni adjustment instead controls family-wise error, the probability of at least one false rejection, and is often more conservative when many hypotheses are tested. Under arbitrary dependence, a more conservative procedure such as Benjamini–Yekutieli can control FDR under broader conditions. Choice of procedure depends on the error criterion and assumptions; no adjustment makes poorly calibrated raw p-values valid.
Finally, statistical evidence and biological importance are distinct. A very small effect can be estimated precisely with a large sample and yield a small adjusted p-value. A potentially large effect can have a wide interval and fail to pass a threshold because the data are imprecise. Neither outcome alone settles the scientific question. Interpret magnitude, direction, uncertainty, design quality, and biological context together, and avoid causal language unless the design and assumptions justify it.
Check Your Understanding
Answer each question, then compare with the solutions below.
- In a family of 100 gene tests, what does a BH adjusted p-value of 0.03 mean, and what does it not mean?
- A gene has log2 fold change \(-1\). What is its estimated fold change, and how should its direction be described?
- For sorted p-values 0.01, 0.04, and 0.20, calculate the BH adjusted p-values.
- Why can a gene with a very small adjusted p-value still have little biological importance?
Solutions. (1) The gene would be rejected by BH at FDR levels at least 0.03 for that specified family and analysis. It is not a 3% probability that the null is true, nor a guarantee that 3% of the realized rejections are false. (2) \(2^{-1}=0.5\): the estimated expression is half as large in the numerator condition as in the reference condition, a 50% decrease. (3) Here \(m=3\). The unmonotonized values are \(3(0.01)/1=0.03\), \(3(0.04)/2=0.06\), and \(3(0.20)/3=0.20\); they are nondecreasing, so the adjusted values are 0.03, 0.06, and 0.20. (4) A large, precise data set can provide strong evidence for a small nonzero effect. The adjusted p-value measures evidence under multiple testing, not the size or practical importance of the effect; inspect the fold change and its uncertainty in context.