Adjusting for Multiple Pathway Tests

Why Adjust Pathway P-Values?

A pathway analysis can produce a p-value for each of many gene sets. Even if every pathway null hypothesis is true, a collection of tests can yield small p-values just by chance. The concern is not resolved by examining each pathway at the same unadjusted significance level: the more tests performed, the more opportunities there are for chance findings. False-discovery-rate adjustment addresses this multiplicity while targeting the expected proportion of false findings among the pathways declared significant.

The p-values must come from valid pathway tests, with their null hypotheses and analysis choices specified. Adjustment does not repair invalid p-values, inappropriate gene-level inputs, or a test family chosen after inspecting the results. The earlier material on competitive and self-contained gene-set tests distinguishes the kinds of null hypotheses a pathway p-value can represent; the adjustment described here applies to either kind, provided the p-values are valid for their respective tests.

Definition: For a prespecified family of $m$ tests, let $R$ be the number of rejected null hypotheses and $V$ the number of rejected null hypotheses that are in fact true. The false discovery rate is $$ \operatorname{FDR}=E\left[\frac{V}{\max(R,1)}\right]. $$ The denominator convention makes the fraction zero when $R=0$. A procedure controls the FDR at level $\alpha$ if this expectation is at most $\alpha$, under its stated assumptions.

FDR control is a repeated-sampling error-rate statement about a procedure applied to a family of tests. It does not say that a particular set of reported pathways has a specified probability of containing a false result. It is also different from family-wise error control, which limits the probability of making even one false rejection.

The Benjamini–Hochberg Procedure

The Benjamini–Hochberg (BH) procedure is a standard FDR method. Begin with $m$ valid p-values and sort them from smallest to largest, writing $p_{(1)}\leq p_{(2)}\leq\cdots\leq p_{(m)}$. For a chosen target level $\alpha$, compare each ordered p-value with its rank-specific threshold $i\alpha/m$. Find the largest rank $k$ for which $p_{(k)}\leq k\alpha/m$. Reject the hypotheses at ranks $1,\ldots,k$. If no rank satisfies the inequality, reject none.

Under independence of the tests, BH controls FDR at or below $\alpha$. It also controls FDR under certain positive dependence conditions, commonly expressed as positive regression dependency on the subset of true null hypotheses (PRDS). Overlapping pathways can create dependence, but overlap alone does not establish PRDS. Thus the procedure can always be calculated, while the formal guarantee depends on the joint behavior of the p-values.

p_{(i)}\leq \frac{i\alpha}{m}, \qquad k=\max\left\{i:p_{(i)}\leq\frac{i\alpha}{m}\right\}.

An equivalent way to report the BH results is with adjusted p-values, often called BH q-values in software output. For the p-value at ordered rank $i$, calculate $m p_{(j)}/j$ for ranks $j\geq i$, take the smallest of those values, and cap the result at 1. This reverse cumulative minimum is essential: simply multiplying every p-value by $m/i$ can produce adjusted values that decrease as the original p-values increase.

q_{(i)}=\min\left(1,\ \min_{j\geq i}\frac{m p_{(j)}}{j}\right).

The adjusted p-value is the smallest BH target level at which that hypothesis would be rejected, given the same family and p-values. At a chosen level $\alpha$, the BH rejection set is equivalently the set with $q_{(i)}\leq\alpha$. The adjusted p-value is not the probability that the pathway null hypothesis is true, nor is it a measure of the size or direction of a pathway effect.

Worked Example: Applying BH at a Five-Percent FDR Level

Suppose six hypothetical pathway tests have p-values $0.001,\ 0.008,\ 0.012,\ 0.030,\ 0.20,\ 0.60$. Assume each p-value is valid, and consider the entire set of six tests as the family. At $\alpha=0.05$, the ordered thresholds are $i(0.05)/6$.

The thresholds for ranks 1 through 6 are $0.00833,\ 0.0167,\ 0.0250,\ 0.0333,\ 0.0417,\ 0.0500$, respectively. The rank-4 p-value satisfies $0.030\leq0.0333$. Rank 5 fails because $0.20>0.0417$, and rank 6 fails because $0.60>0.0500$. Rank 4 is therefore the largest qualifying rank, so BH rejects the first four pathway null hypotheses.

The adjusted values are calculated from the largest ranks backward. At rank 6, $6(0.60)/6=0.60$, so $q_{(6)}=0.60$. At rank 5, the scaled value is $6(0.20)/5=0.24$, and the cumulative minimum is $q_{(5)}=0.24$. At rank 4, $6(0.030)/4=0.045$, giving $q_{(4)}=0.045$. At rank 3, $6(0.012)/3=0.024$, so $q_{(3)}=0.024$. At rank 2, $6(0.008)/2=0.024$, and at rank 1, $6(0.001)/1=0.006$. The adjusted values in rank order are therefore $0.006,\ 0.024,\ 0.024,\ 0.045,\ 0.24,\ 0.60$. The first four are at most 0.05, matching the step-up decision.

The conclusion is that these four tests meet the BH rejection rule at a five-percent target, subject to the procedure's assumptions. It is not a claim that each finding has a five-percent chance of being false.

Why the Adjusted-P-Value Formula Works

The step-up rule and the adjusted-p-value rule give the same decisions. To see why, fix a rank $i$. Its adjusted value is at most $\alpha$ exactly when there is some rank $j\geq i$ such that $m p_{(j)}/j\leq\alpha$. Rearranging gives $p_{(j)}\leq j\alpha/m$, meaning that rank $j$ qualifies in the step-up comparison. If such a $j$ exists, the largest qualifying rank $k$ is at least $i$, and the BH procedure rejects rank $i$. Conversely, if rank $i$ is rejected, the largest qualifying rank $k$ satisfies $k\geq i$; using $j=k$ in the minimum shows $q_{(i)}\leq m p_{(k)}/k\leq\alpha$. Thus the two methods select exactly the same ranks.

Worked Example: Interpreting an Adjusted Value

Consider five hypothetical pathway tests with ordered p-values $0.004,\ 0.015,\ 0.020,\ 0.21,\ 0.40$. For rank 3, the scaled values from ranks 3 through 5 are $5(0.020)/3=0.0333$, $5(0.21)/4=0.2625$, and $5(0.40)/5=0.40$. Their minimum is $q_{(3)}=0.0333$.

This pathway would be rejected by BH at $\alpha=0.05$, holding the five-test family and its p-values fixed. The interpretation is that its smallest BH rejection target is about 0.0333. It does not mean the pathway's null has a 3.33% probability of being true, nor does it say how large the pathway-level association is. The pathway's test statistic, estimated effect or enrichment magnitude, and biological interpretation remain separate evidence to report.

Overlapping Pathways and Dependence

Pathways often share genes. If two pathway statistics are influenced by many of the same genes, their p-values may be dependent: a change in the shared genes can move both tests. Related pathway definitions, common sample-level variation, and the structure of the gene-level statistics can also contribute to dependence. The earlier material on overlapping pathways describes why significant results from related gene sets are not necessarily independent biological discoveries.

Overlap does not change the BH arithmetic: sort the p-values, use the total number $m$ in the family, and compute the same thresholds and adjusted values. It does matter for how strongly one can rely on the formal FDR guarantee. One must not infer that BH's independence condition holds simply because the p-values are listed separately, or that all overlap produces dependence in the same direction. When the dependence structure does not meet a known BH condition, the usual BH guarantee is not automatic.

Worked Example: A Dependence-Robust Alternative

Return to the six p-values $0.001,\ 0.008,\ 0.012,\ 0.030,\ 0.20,\ 0.60$, and suppose the pathways share genes in a way that makes the dependence structure difficult to justify as satisfying BH's positive-dependence condition. A conservative alternative for arbitrary dependence is the Benjamini–Yekutieli (BY) procedure. It multiplies the BH thresholds by $1/H_m$, where $H_m=\sum_{j=1}^m1/j$, or equivalently multiplies BH adjusted p-values by $H_m$, capping at 1.

For $m=6$, $H_6=1+1/2+1/3+1/4+1/5+1/6=2.45$. The BH adjusted value for the rank-4 pathway was $0.045$, so its BY adjusted value is $2.45(0.045)=0.1103$. The BH values $0.006,\ 0.024,\ 0.024,\ 0.045,\ 0.24,\ 0.60$ become $0.0147,\ 0.0588,\ 0.0588,\ 0.1103,\ 0.588,\ 1.00$, after capping at 1. At a 0.05 target, only the rank-1 pathway is rejected by BY.

This difference reflects a more conservative procedure, not a change in the observed pathway p-values or in the biological strength of the evidence. BY provides an FDR guarantee under arbitrary dependence for valid p-values, whereas ordinary BH is more powerful when its dependence conditions are credible. The choice should be justified by the intended error control and dependence assumptions rather than by which method produces more findings.

Choosing and Reporting the Test Family

The family is the collection of hypotheses to which the adjustment is applied. Its definition affects every adjusted value: a p-value that is small within a set of 20 tests can receive a different adjustment when included in a family of 2,000. Define the family in relation to the scientific question and the analysis plan, not after selecting the most attractive results. For example, if several pathway databases, gene-set sizes, or analysis variants were searched and only one set of results is reported, adjusting only the reported subset can understate the multiplicity involved in the search.

A clear report states the test family, number $m$ of tests, adjustment procedure, target level, and the adjusted results. It should also identify the pathway test and null hypothesis, describe relevant gene-set overlap or dependence concerns, and report evidence about magnitude or direction separately when available. An adjusted p-value is a multiplicity-aware measure of incompatibility with a null hypothesis under the specified testing and adjustment procedure; it is not an effect size, a measure of biological importance, or proof of causality.

Common Mistakes and Diagnostics

Using the wrong denominator. The BH rank threshold uses the number $m$ of tests in the stated family, not the rank, number of selected pathways, or number of pathways with small p-values.
Stopping at the first failed threshold. BH is a step-up procedure. Find the largest qualifying rank; a later p-value can qualify even if an earlier rank does not.
Reporting raw scaled values as adjusted p-values. The adjusted values require the reverse cumulative minimum and a cap at 1, so they are nondecreasing with ordered p-values.
Calling q a probability that a result is false. It is an adjusted rejection threshold under the procedure, not a posterior probability for an individual pathway.
Treating overlapping significant pathways as independent confirmations. Shared genes can yield related results. Report overlap and interpret related findings as potentially connected, not as independent replications.
Assuming adjustment validates the analysis. Check that the original pathway p-values are valid, that the test family is appropriate, and that dependence assumptions for the chosen procedure are defensible.

Useful checks include verifying that adjusted values lie between 0 and 1 and are nondecreasing in p-value rank, confirming that the reported rejection set agrees with $q\leq\alpha$, and checking that all intended tests were included. If those checks fail, common causes include adjusting a subset unintentionally, sorting or mapping p-values incorrectly, or omitting the reverse cumulative minimum.

Key takeaway: BH adjustment controls the expected false-discovery proportion under independence and certain positive dependence conditions. It accounts for the number of pathway tests, but does not measure pathway effect size, establish biological importance, or make overlapping findings independent.

Check Your Understanding

Answer the questions, then compare your reasoning with the solutions.

There are 10 valid pathway p-values. At $\alpha=0.05$, what is the BH threshold for rank 4?
In a family of 8 tests, a pathway has BH adjusted p-value $0.041$. What does this say at target level $0.05$, and what does it not say about the pathway null?
Why can overlapping gene sets complicate the guarantee for ordinary BH adjustment?
What must be done to the scaled BH values $m p_{(i)}/i$ to obtain adjusted p-values?
What is one reason to choose BY rather than BH?

Solutions. (1) The threshold is $4(0.05)/10=0.020$. Rank 4 qualifies only if its ordered p-value is at most 0.020. (2) It is rejected by BH at 0.05 for the specified family; the value 0.041 is not the probability that its null is true or that this individual finding is false. (3) Shared genes can make test statistics and p-values dependent, and overlap alone does not establish the positive-dependence condition under which BH controls FDR. (4) Starting from the largest rank, take reverse cumulative minima of the scaled values and cap at 1; then map the results back to their pathways. (5) BY may be chosen when arbitrary dependence must be accommodated for valid p-values and the assumptions for BH's formal guarantee cannot be justified, accepting that BY is generally more conservative.