# Multiple Testing Correction in Enrichment Analysis

The one-sentence difference: a family-wise correction such as Bonferroni protects you from even a single false positive across every gene set you tested, while a false discovery rate correction such as Benjamini-Hochberg allows a small, quantified fraction of false positives among the gene sets you call significant. Rule of thumb: use Benjamini-Hochberg false discovery rate control for routine over-representation analysis across hundreds or thousands of gene sets, reserve family-wise corrections for small pre-specified pathway panels or confirmatory claims, and use permutation-based false discovery rates when you run gene set enrichment analysis on a full ranked list.

This article explains every correction option you will meet in enrichment analysis, compares them on the criteria that actually decide a choice, works through real numeric examples and published case studies, and gives you a reproducible R workflow plus a troubleshooting list for the failures that come up most often.

## Why Enrichment Analysis Creates a Multiple Testing Problem

Enrichment analysis asks one question many times. An over-representation analysis (ORA) takes a list of differentially expressed genes, tests each gene set in a database such as Gene Ontology or KEGG for an unusually high number of hits, and returns one p-value per gene set. A typical human [Gene Ontology Biological](/knowledge/molecular-biology/gene-ontology-biological) Process run tests several thousand terms in a single function call. If you test 5,000 terms and flag everything with a raw p-value below 0.05, you expect roughly 250 of those flags to be pure noise even if your data contained no biology at all.

Multiple testing correction fixes this by re-scoring each raw p-value in light of how many tests you ran. The corrected number is called an adjusted p-value, an FDR-adjusted p-value, or a q-value depending on the method, and it is the number you report and threshold.

Two recent studies show why this matters in practice. In an analysis of vascular stress signatures before VEGF inhibitor treatment, no individual gene remained significant after multiple testing correction, so the authors placed their emphasis on coordinated biological responses detected through Gene Set Enrichment Analysis (GSEA), which revealed oxidative stress, inflammatory response, and NF-κB signaling enrichment at the pathway level [1]. This is a textbook example of the correction burden: with tens of thousands of genes tested individually, single-gene significance becomes nearly impossible, while the aggregate pathway signal stays detectable. The same pattern appears in [comparative genomics](/blog/guides/comparative-genomics), where three desiccation-tolerance orthogroups showed nominally greater counts of trend-conforming paralogs with uncorrected p below 0.05, and none of the three survived multiple testing correction, a finding the authors interpreted as biologically meaningful but statistically unconfirmed [2].

Both examples carry the same lesson. Enrichment analysis exists partly because the multiple testing burden at the single-feature level is so severe, but enrichment analysis creates its own burden across gene sets, and that burden needs its own correction.

## The Vocabulary That Decides Everything

Before comparing methods, you need six terms defined precisely, because the arguments between correction methods are arguments about which of these quantities they control.

**The family** is the collection of tests you correct across. In enrichment analysis the family is the set of gene sets you tested, not the set of genes. If your ORA tested 7,000 GO terms, your family has 7,000 members. Defining the family correctly is the single most consequential decision in the whole procedure, and getting it wrong is the most common analytical mistake in published enrichment analyses.

**The family-wise error rate (FWER)** is the probability of making even one false positive call across the entire family. Bonferroni and Holm control FWER. At a threshold of 0.05, FWER control means a 5 percent chance that any result in your table is a false positive.

**The false discovery rate (FDR)** is the expected proportion of false positives among the results you declare significant. Benjamini-Hochberg, Benjamini-Yekutieli, Storey q-values, and permutation FDR control FDR. At a threshold of 0.05, FDR control means that if you call 40 pathways significant, you expect about 2 of them to be false.

**A raw p-value** is the ordinary p-value from a single test, computed without reference to the other tests.

**An adjusted p-value** is the raw p-value after correction, reported so that thresholding it at your chosen alpha delivers the stated error guarantee.

**A q-value** is either a synonym for an FDR-adjusted p-value (as used in the GSEA software and the fgsea package) or, in the stricter Storey sense, an estimate of the minimum FDR at which this specific test would count as significant.

One more term belongs in this list because it silently changes results. **The universe (or background)** is the full set of genes that could have been detected in your experiment. An ORA p-value is computed as a contingency test against this universe, so using all annotated genes when your platform only measures 20,000 of them inflates significance across every gene set at once.

## Multiple Testing Correction Options at a Glance

The table below is the core comparison. Every enrichment tool you will touch offers some subset of these options, and the differences between them are the differences between publishable and unpublishable conclusions.

| Method | Error rate controlled | How it adjusts | Dependence assumption | Power | Typical role in enrichment analysis |
|--|--|--|--|--|--|
| Bonferroni | FWER | Multiplies each raw p by the number of tests | None (valid under any dependence) | Lowest | Small pre-specified panels, extra-conservative confirmation |
| Holm | FWER | Step-down comparison of sorted p-values to shrinking thresholds | None | Higher than Bonferroni at the same FWER | Confirmatory follow-ups on a short pathway list |
| Benjamini-Hochberg (BH) | FDR | Step-up comparison of sorted p-values to scaled thresholds | Independence or positive regression dependence | High | Default for ORA in clusterProfiler, g:Profiler, and most R pipelines |
| Benjamini-Yekutieli (BY) | FDR | BH with a harmonic-number penalty for dependence | Arbitrary dependence | Lower than BH | Pathway collections with heavy overlap or strong correlation |
| Storey q-value | FDR | Estimates the proportion of true nulls and reweights BH | Similar to BH, with tuning | Higher than BH when many true signals exist | Large GO or KEGG screens where many pathways are genuinely affected |
| Permutation-based (GSEA phenotype or gene-set permutation) | Empirical FDR or FWER | Builds an empirical null by resampling class labels or gene sets | Exchangeable samples | Depends on sample size | Ranked-list GSEA where no per-gene p-values exist |
| Covariate-weighted FDR (for example, Independent Hypothesis Weighting) | FDR | Assigns weights using an informative covariate such as set size | Covariate independent of p under the null | Higher than unweighted BH | Very large screens with an informative auxiliary variable |

## The Correction Methods Explained

Each option below gets its own definition, its pros and cons, and its honest scope of use.

### Bonferroni Correction

The Bonferroni correction multiplies every raw p-value by the number of tests (or, equivalently, declares significance when the raw p falls below alpha divided by the number of tests). It guarantees that the probability of even one false positive across the family stays at or below your alpha level, and it requires no assumption about how the tests correlate.

Pros. It is simple, exact, and bulletproof under any dependence structure. When a result survives Bonferroni in an enrichment analysis, almost nobody challenges it.

Cons. It is brutally conservative when the family is large. With 7,000 GO terms, a raw p of 0.00001 becomes 0.07 and misses the cutoff, which throws away real signal. It also ignores a statistical fact that helps you: pathway p-values are positively correlated because gene sets overlap, so treating them as fully independent extremes wastes power.

Bonferroni fits enrichment analysis in narrow but real situations. A study of PSEN1 missense variant pathogenicity used Fisher's exact tests with Bonferroni correction to assess enrichment of functional features across a small set of defined comparisons, a design where the family is small and the claim is confirmatory [3]. Use Bonferroni when you test tens of gene sets, not thousands, or when your journal or reviewer explicitly asks for family-wise control.

### Holm Correction

The Holm (step-down Bonferroni) procedure sorts raw p-values from smallest to largest and compares the smallest to alpha divided by m, the second smallest to alpha divided by (m - 1), and so on, rejecting until the first failure. The adjusted p-value form of Holm is what p.adjust(method = "holm") returns.

Pros. It controls FWER under any dependence, exactly like Bonferroni, and it is uniformly more powerful, sometimes meaningfully so on small pathway panels.

Cons. It is still a family-wise method, so on large GO screens it remains far too conservative for exploratory work. In practice, Holm appears in enrichment analysis less often than either Bonferroni or BH, usually because analysts who want FWER control default to Bonferroni for its simplicity and analysts who want power default to BH.

### Benjamini-Hochberg (BH)

The BH procedure sorts the raw p-values from smallest to largest, compares each to the threshold (i / m) × alpha for its rank i, finds the largest rank that passes, and declares all tests at or above that rank significant. Its adjusted p-values have the convenient property that thresholding them at 0.05 gives exactly the BH calls at alpha of 0.05.

BH controls the false discovery rate under independence and, importantly, under positive regression dependence, which is the situation pathway p-values usually occupy because overlapping gene sets induce positive correlation. This is why BH became the default in essentially every enrichment package. clusterProfiler's enrichGO and enrichKEGG functions use BH as their default pAdjustMethod, and fgsea returns BH-adjusted padj values alongside raw p-values.

Pros. High power. Interpretation is intuitive and quantitative: at FDR 0.05 you expect about 5 percent of your called pathways to be false. It matches how omics results are actually used, where a ranked list of 30 plausible pathways is more useful than a list of 2 ultra-safe ones.

Cons. Under the null with no true signal at all, calling anything significant means those calls are all false, and a fraction of 0.05 of a small call list can still be one misleading pathway. BH also assumes positive dependence, and adversarial or exotic dependence structures fall outside its guarantee, which is where Benjamini-Yekutieli enters.

Published enrichment analyses use BH constantly. A prospective study of serum microRNAs in relapsing-remitting multiple sclerosis identified eight differentially expressed miRNAs between relapse and remission after Benjamini-Hochberg correction, with hsa-miR-122-5p showing the largest fold change (fold change = 3.0, adjusted p = 0.043) [4]. A lung tissue study associated fourteen miRNA loci with quantitative emphysema at BH-FDR below 0.05, using linear regression with heteroskedasticity-robust standard errors as the underlying test [5]. And a rare-variant analysis of severe COVID-19 reported five genes associated with severe disease at FDR below 0.05 via the variable threshold method, after classical single-variant and burden tests found nothing following correction [6]. Note the range of contexts there: BH serves qPCR panels, microRNA linear models, and gene-level burden tests alike, because the correction operates on whatever p-values the family produces.

### Benjamini-Yekutieli (BY)

BY is BH with an extra penalty factor, the sum of the harmonic series up to m (approximately ln(m) + 0.577 for large m), which buys validity under arbitrary dependence between tests. At m = 5,000 the penalty factor is roughly 9.1, so every BH threshold shrinks by about ninefold.

Pros. The FDR guarantee holds no matter how correlated the tests are, which makes it defensible for pathway collections where overlap is extreme or where a critic questions the positive-dependence assumption.

Cons. The power loss is severe on large families, and in most omics settings it is more conservative than the data require, because pathway test statistics generated from the same gene measurements are overwhelmingly positively dependent. Use BY when you specifically need dependence-free FDR control, for example when reviewers question BH validity for a small, heavily overlapping curated panel, not as a default.

### Storey q-Value

Storey's method estimates π₀, the proportion of tests whose null hypothesis is true, then uses that estimate to convert raw p-values into q-values with more power than BH when many tests are truly non-null. The qvalue package in Bioconductor implements it, and the estimate depends on a tuning parameter (lambda) that you should inspect rather than accept blindly.

Pros. When a screen genuinely hits dozens or hundreds of pathways, π₀ is well below 1, and Storey q-values shrink relative to BH, recovering real hits that BH misses.

Cons. The π₀ estimate is unstable when true signals are sparse or when the p-value histogram is skewed by correlated tests, which is common in pathway panels. An over-optimistic π₀ inflates false discoveries exactly when you can least detect it. Storey q-values work best on large, well-behaved p-value distributions, and they require you to report the π₀ estimate so readers can judge the calibration.

### Permutation-Based Correction in GSEA

Gene Set Enrichment Analysis does not start from per-gene p-values. It walks down a ranked list of all genes (ranked by a statistic such as the signal-to-noise ratio or the DESeq2 stat) and computes an enrichment score for each gene set, then derives significance by permuting. Phenotype permutation shuffles the class labels and recomputes the whole ranked list each time, preserving gene-gene correlation. Gene-set permutation shuffles gene labels instead, which is faster but destroys that correlation structure. The GSEA documentation recommends phenotype permutation when you have roughly seven or more samples per class and falls back to gene-set permutation when you do not.

After normalizing enrichment scores across gene sets and permutations, GSEA reports a nominal p-value and an FDR q-value estimated from the permutation null distribution across all gene sets and permutations. This is a resampling-based FDR estimate rather than a BH adjustment, and it is the correction framework that pathway-level enrichment studies typically report.

Pros. It accounts for gene-gene correlation directly instead of assuming dependence structure, and it uses the whole ranked list, which captures weak coordinated shifts that ORA on a cut gene list cannot see. The VEGF inhibitor vascular study worked exactly this way, using GSEA on pre-treatment tumor [RNA sequencing](/blog/guides/rna-sequencing) data to identify pathway-level signatures after individual genes failed multiple testing correction [1].

Cons. Phenotype permutation needs adequate sample sizes, and with small cohorts the empirical null is coarse, so q-values are unstable and small q-values are unattainable. Gene-set permutation, the small-sample fallback, treats genes as exchangeable and can be anti-conservative for correlated sets. Permutation is also slower and its results are stochastic unless you fix the random seed, which matters for reproducibility, a point that integrated frameworks for reproducible enrichment workflows now treat as a first-class concern [7].

### Covariate-Weighted FDR

Independent Hypothesis Weighting (IHW) assigns each test a weight derived from a covariate that is informative about power but independent of the p-value under the null (for ORA, a natural covariate is gene set size or the number of genes in the intersect). Weighted BH then controls FDR at the nominal level while spending more of the error budget on tests likely to be true discoveries. The ihw function in the IHW package fits the weights.

Pros and cons in brief. IHW raises power on large screens and is statistically principled, but it adds a layer of machinery that most enrichment workflows skip, and a poorly chosen covariate can break the null-independence requirement. Treat it as an advanced option rather than a default, and mention in your methods that weighting was applied and which covariate was used.

## Worked Example: The Same Ten Pathways Under Three Corrections

![Formula panel showing Bonferroni p times m and Holm p times m minus rank plus one, worked for raw p 0.0004 at rank 2 of 10.](/article-images/illustrations/multiple-testing-correction-enrichment-analysis-bonferroni-holm-example.webp)

Nothing communicates the trade-offs like the same p-values run through three corrections. Suppose an ORA tested 10 gene sets (a small curated panel, kept small so the arithmetic stays visible) and produced these raw p-values, already sorted from smallest to largest.

| Rank | Raw p | Bonferroni (p × 10) | Holm adjusted | BH adjusted | Passes BH at 0.05 |
|--|--|--|--|--|--|
| 1 | 0.0001 | 0.001 | 0.001 | 0.001 | Yes |
| 2 | 0.0004 | 0.004 | 0.0036 | 0.002 | Yes |
| 3 | 0.0009 | 0.009 | 0.0072 | 0.003 | Yes |
| 4 | 0.0020 | 0.020 | 0.014 | 0.005 | Yes |
| 5 | 0.0048 | 0.048 | 0.0288 | 0.0096 | Yes |
| 6 | 0.0100 | 0.100 | 0.050 | 0.0167 | Yes |
| 7 | 0.0200 | 0.200 | 0.080 | 0.0286 | Yes |
| 8 | 0.0400 | 0.400 | 0.160 | 0.050 | Yes |
| 9 | 0.0600 | 0.600 | 0.240 | 0.0667 | No |
| 10 | 0.0900 | 0.900 | 0.360 | 0.09 | No |

The conclusions differ sharply. Bonferroni and Holm pass the top five gene sets at alpha 0.05. BH passes the top eight, including rank 8 exactly at 0.050. The Holm adjusted values in this table are computed as p(rank i) × (10 - i + 1) with monotonicity enforced, which is why rank 6 sits exactly at 0.05 and passes while rank 7 fails. If this panel fed a confirmatory claim, the family-wise methods would be the defensible choice and you would report five hits. If this were an exploratory screen feeding follow-up experiments, BH would be the defensible choice and you would carry eight candidates forward, knowing roughly 5 percent of them (fewer than one, in expectation) might be false.

You can reproduce every number in this table with one line of base R.

```r
p <- c(0.0001, 0.0004, 0.0009, 0.002, 0.0048, 0.01, 0.02, 0.04, 0.06, 0.09)
p.adjust(p, method = "bonferroni")
p.adjust(p, method = "holm")
p.adjust(p, method = "BH")
p.adjust(p, method = "BY")
```

The BY output at this family size is nearly identical to BH because the harmonic penalty is small when m is small. Rerun the same p-values with 5,000 appended weak signals and the BY penalty becomes visible, which is exactly how the method behaves on real GO screens.

## Worked Scenarios From Real Studies

### Scenario 1: No Genes Survive, Pathways Do

A transcriptomic study of vascular toxicity from VEGF inhibitors analyzed pre-treatment tumor RNA sequencing data with GSEA, and no individual gene remained significant after multiple testing correction [1]. Rather than treating this as a failure, the analysis shifted its unit of inference to coordinated biological responses and found significant baseline enrichment of oxidative stress, inflammatory response, and NF-κB signaling pathways in patients who later received VEGF inhibitors, which functional validation in human vascular cells then supported through reactive oxygen species production, cytokine release, and extracellular vesicle release [1]. The methodological point: when the single-feature family is huge and effects are modest and distributed, pathway-level inference with a permutation FDR is often the only design with adequate power, not a consolation prize.

### Scenario 2: Two-Tier Thresholding in Metabolomics

An untargeted metabolomics comparison of Marselan and Cabernet Sauvignon grape berries first identified 236 candidate differential metabolites using VIP above 1 and unadjusted p below 0.05, then applied a stricter tier requiring FDR at or below 0.05 and an absolute log2 fold change of at least 1, which retained 34 features, from which 13 biologically interpretable features were selected for presentation after review of annotation evidence [8]. In the downstream KEGG enrichment, flavonoid biosynthesis remained significant after multiple testing correction, while phenylalanine metabolism and flavone and flavonol biosynthesis showed distinct differential metabolite mapping patterns [8]. The methodological point: layered thresholds (raw screen, then FDR, then effect size and interpretability) are a legitimate and transparent workflow, provided each tier is stated in the methods and the uncorrected tier is labeled as a screen.

### Scenario 3: BH in Small Molecular Panels

The multiple sclerosis serum miRNA study ran a targeted qPCR panel with paired statistical testing and Benjamini-Hochberg correction, and reported an adjusted p of 0.043 for its strongest miRNA [4]. With a small targeted family, BH retains enough power to reach conventional significance while still carrying an FDR interpretation. The methodological point: panel studies are the sweet spot where BH is both powerful and unimpeachable, and where a family-wise method would often have erased the finding entirely.

### Scenario 4: Honest Null Results After Correction

Three studies in the recent literature reported that their primary comparisons produced nothing significant after correction and said so plainly. The CHMP4C colorectal cancer study found that neither MMR nor MSI status was significantly associated with CHMP4C alterations after multiple-testing correction, even though uncorrected counts showed directional representation [9]. The rare-variant COVID study reported that classical single-variant and burden tests identified no significant associations after multiple testing correction before pivoting to network-based clustering, and separately reported five genes at FDR below 0.05 from the variable threshold analysis [6]. And the resurrection plant comparative genomics study reported that its three nominally significant orthogroups (dehydrin, a NAC-domain transcription factor, and a LEA hydroxyproline-rich glycoprotein) failed correction while arguing their biological plausibility explicitly as a hypothesis for future work [2]. The methodological point: a null after correction is a result. Reporting it honestly, with the uncorrected values visible and labeled, is better science than silently switching to nominal p-values.

## A Reproducible R Workflow

The following snippets cover the three correction contexts you will encounter most: manual adjustment, ORA, and ranked-list GSEA. Each is runnable as written, assuming standard Bioconductor packages.

First, the universal tool. p.adjust ships with base R and implements bonferroni, holm, hochberg, hommel, BH, BY, and fdr (which is an alias for BH).

```r
raw_p <- c(0.001, 0.012, 0.049, 0.20, 0.63)
p.adjust(raw_p, method = "BH")
data.frame(raw = raw_p,
           bonferroni = p.adjust(raw_p, "bonferroni"),
           holm = p.adjust(raw_p, "holm"),
           BH = p.adjust(raw_p, "BH"),
           BY = p.adjust(raw_p, "BY"))
```

Second, over-representation analysis with clusterProfiler. The default pAdjustMethod is BH, and the default cutoffs are strict, which surprises people whose results vanish.

```r
library(clusterProfiler)
library(org.Hs.eg.db)

ego <- enrichGO(
  gene          = sig_genes,
  OrgDb         = org.Hs.eg.db,
  keyType       = "ENTREZID",
  ont           = "BP",
  pAdjustMethod = "BH",
  pvalueCutoff  = 0.05,
  qvalueCutoff  = 0.2,
  readable      = TRUE
)

as.data.frame(ego)
```

A practical debugging trick: set both cutoffs permissively (pvalueCutoff = 1, qvalueCutoff = 1), inspect the full table with raw and adjusted columns side by side, then filter in a separate step. This separates "the enrichment test found nothing" from "the correction erased everything," which are different problems with different fixes.

Third, ranked-list GSEA with fgsea, which computes permutation-based enrichment scores and reports BH-adjusted padj values.

```r
library(fgsea)

genes_ranked <- sort(de_stats$stat, decreasing = TRUE)
names(genes_ranked) <- de_stats$entrez

fgs <- fgsea(pathways = msigdb_pathways,
             stats    = genes_ranked,
             minSize  = 15,
             maxSize  = 500)

fgs <- fgs[order(pval)]
fgs[, .(pathway, NES, pval, padj, size)]
```

Finally, for Storey q-values on a large screen.

```r
library(qvalue)
qobj <- qvalue(p = ora_pvalues, lambda = seq(0.05, 0.95, 0.05))
qobj$pi0
summary(qobj)
ora_qvalues <- qobj$qvalues
```

Report the π₀ estimate whenever you use Storey q-values, because its value determines how much the results lean on the tuning.

## How to Choose the Best Correction for Your Design

![Decision path asking about family size, permutation null and claim strength, ending in BH, Bonferroni or Holm, or permutation FDR.](/article-images/illustrations/multiple-testing-correction-enrichment-analysis-correction-choice.webp)

There is no single best multiple testing correction in enrichment analysis. There is a best correction for a given family size, dependence structure, and claim strength. Work through the following list in order.

1. How many gene sets are in the family? Under about 50, family-wise control is affordable and Bonferroni or Holm are defensible for confirmatory claims. Above about 200, family-wise control discards most real signal and FDR is the appropriate framework.
2. Are the gene sets heavily overlapping? GO Biological Process branches overlap massively, which induces positive dependence. BH handles positive dependence, so overlap is not a reason to abandon BH, but it is a reason to avoid pretending the tests are independent.
3. Does your test produce its own empirical null? If you run phenotype-permuted GSEA with adequate sample sizes, use its built-in FDR q-value. Do not additionally apply BH on top of permutation q-values.
4. What is the claim? Exploratory hypothesis generation licenses FDR control at 0.05 or even 0.10 with explicit labeling. A claim that a specific pathway is dysregulated as a headline result deserves family-wise or replicated support.
5. What will reviewers in your field expect? In human genomics and transcriptomics, BH at FDR 0.05 is the lingua franca, and deviations need justification. In hypothesis-driven panels, Bonferroni appears routinely, as in the PSEN1 variant enrichment study [3].

Applied to real designs: for a standard bulk RNA-seq ORA with thousands of GO terms, use BH. For a GSEA on a ranked list with 20 or more samples per class, use phenotype permutation with the GSEA or fgsea FDR. For a 40-miRNA qPCR panel, use BH, as the multiple sclerosis study did [4]. For a GWAS-derived pathway or cell-type enrichment with very precise p-values, use FDR thresholds and label single-instrument or underpowered estimates as exploratory regardless of their FDR, exactly as the atopic dermatitis and asthma cross-trait analysis did when it called single-instrument estimates exploratory no matter the false discovery rate [10]. For Mendelian randomization screens across many exposures, FDR correction is standard, as in the serum metabolite and tinnitus analysis that applied FDR correction across 486 metabolites [11].

## Common Mistakes and Limitations

This section doubles as a bioinformatics troubleshooting checklist for multiple testing correction in enrichment analysis, covering the failures that produce wrong results, silent inconsistencies, and reviewer objections.

1. **Correcting the wrong family.** The most frequent conceptual error is correcting across genes before running enrichment, or believing that per-gene FDR carries over to the pathway table. The p-values you report in the enrichment table are the ones that need the correction, and the family is the set of gene sets tested. A related error is forgetting that DESeq2 and edgeR already apply some internal multiplicity thinking (DESeq2 performs independent filtering to increase power), so the gene-level padj you started from is itself a corrected quantity.

2. **Applying correction twice.** Some pipelines adjust gene p-values for FDR, run ORA on the survivors, then adjust ORA p-values again with a stricter threshold "to be safe." Double correction is not extra-safe, it distorts the error rate in an uncontrolled direction and destroys power. Choose one correction per reported table.

3. **Background (universe) mismatch.** If your RNA-seq pipeline only quantifies genes expressed above a floor, but the ORA compares your list against all annotated genes, every pathway containing lowly expressed genes gets a biased test. Use the set of genes that could have been detected as the universe (clusterProfiler accepts a universe argument), and check that the universe is not silently reset when you switch tools.

4. **Comparing tools without reading their defaults.** clusterProfiler defaults to BH with a qvalueCutoff of 0.2. Other tools default differently. Two pipelines given identical input can produce different pathway lists purely from default cutoffs and universe handling, and a benchmarking framework built to standardize differential analysis, pathway enrichment, and visualization across heterogeneous tools exists precisely because this fragmentation compromises reproducibility and interpretability [7]. When results differ between tools, enumerate defaults before suspecting the data.

5. **Permutation with too few samples.** Phenotype permutation in GSEA needs enough samples to build a usable empirical null. With small cohorts, gene-set permutation is the fallback, and its weakness against gene-gene correlation should temper the strength of your conclusions. If your cohort is small, treat permutation q-values near 0.05 as provisional and seek replication.

6. **Propagating upstream RNA-seq errors into the family.** RNA-seq errors upstream of enrichment, including batch effects, inconsistent normalization, poor-quality samples, and mislabeled or duplicated identifiers, change the input gene list and therefore change which tests exist and how correlated they are. A differential expression step that ignores batch structure can produce a gene list whose pathway hits are artifacts of the batch, no matter how carefully you correct the pathway p-values afterward. Audit QC metrics, normalization choices, and identifier mapping before interpreting any adjusted p-value, and note that gene-gene interaction structure and coexpression bias in differential gene selection can distort downstream enrichment in ways that marginal per-gene models do not capture [7].

7. **Selective reporting across tiers.** Two-tier thresholds are legitimate when declared, but reporting only the tier that produced the result you liked (VIP-based hits without the FDR tier, or nominal pathway p-values when corrected ones failed) is selective inference. State every tier, including the ones that failed, and label exploratory analyses as such.

8. **Interpreting FDR at the level of a single pathway.** An adjusted p of 0.03 for one pathway does not mean that pathway has a 3 percent chance of being false. It means that among the pathways you call significant at this threshold, the expected false proportion is 3 percent. Single-result language ("this pathway has a 97 percent chance of being real") is a misreading that reviewers catch.

9. **Machine-learning feature importance is not enrichment inference.** As a caution from adjacent methodology, a graph attention network re-evaluation on peripheral blood RNA-seq found that attention-derived explanations were not statistically distinguishable from a rigorously matched permutation null even when classification accuracy looked acceptable, and the analysis survived scrutiny because it included permutation testing and orthogonal differential-expression and pathway checks from the start [12]. If your pipeline combines ML feature selection with downstream enrichment, the correction applies to the enrichment tests, and the ML component needs its own null comparison.

10. **Forgetting that correction thresholds are conventions.** FDR 0.05 is a convention, not a law. In exploratory screens with strong downstream validation (functional assays, independent cohorts), 0.10 with explicit labeling can be reasonable. In confirmatory single-cohort claims, 0.05 may still be too loose. Match the threshold to the cost of a false lead in your specific project.

One limitation deserves its own paragraph. Multiple testing correction changes what your results mean, not just how small they look. After BH, your table answers the question "which gene sets would I expect to keep if I repeated this experiment," and after Bonferroni it answers "which gene sets would survive even if I demanded zero false positives overall." Neither question is "which pathways are true," and no correction rescues a design with the wrong universe, unmodeled batch effects, or a gene list contaminated by technical artifacts.

## Frequently Asked Questions

### Should I correct p-values across genes or across pathways in enrichment analysis?

Correct across the tests you are reporting, which in an enrichment table means across gene sets. Per-gene correction happens in the differential expression step (for example, DESeq2 or edgeR adjusted p-values), and the pathway-level family then needs its own correction across the tested gene sets. Applying the pathway-level correction to gene p-values, or assuming gene-level correction covers pathways, is wrong in both directions.

### Is an FDR of 0.05 the same as a p-value of 0.05?

No. A raw p-value of 0.05 means a 5 percent chance of a result this extreme if that one gene set were unaffected. An FDR-adjusted value of 0.05 means that among all gene sets you call significant at this threshold, the expected proportion that are false is 5 percent. The two quantities answer different questions and are not interchangeable, which is why adjusted values are almost always larger than raw ones.

### What is the best multiple testing correction for enrichment analysis?

For routine over-representation analysis across hundreds or thousands of gene sets, Benjamini-Hochberg FDR is the best default and the field standard. For GSEA on a ranked list with adequate sample sizes, use the permutation-based FDR the method provides. For small pre-specified panels or confirmatory claims, Bonferroni or Holm family-wise control is the stronger and more defensible choice.

### Why do clusterProfiler and other tools give me different adjusted p-values for the same gene list?

Differences usually come from the correction method and cutoff defaults, the universe (background) each tool uses, gene identifier mapping, and whether the tool tests over-representation or ranks. Enumerate each tool's defaults, force a common universe and common identifier type, and rerun before concluding that one tool is wrong.

### Can I report raw p-values alongside adjusted p-values in my enrichment table?

Yes, and you should. Most published tables show both, which lets readers see how much of a hit depends on the correction and makes tiered screening transparent. The rule is that significance claims and colored figures key off the adjusted value, with raw values shown for context and labeled clearly.

### My adjusted p-values come back as exactly 1.0. What does that mean?

In Benjamini-Hochberg, adjusted p-values are capped at 1 because a probability cannot exceed 1, so weak p-values far down the ranking get floored at 1. An adjusted value of exactly 1.0 means the test carried essentially no evidence after accounting for the number of tests, not that the calculation failed. If a result you care about shows this, look at its raw p-value and the family size to understand the gap.

<script type="application/ld+json">
{
  "@context": "https://schema.org",
  "@type": "FAQPage",
  "mainEntity": [
    {
      "@type": "Question",
      "name": "Should I correct p-values across genes or across pathways in enrichment analysis?",
      "acceptedAnswer": {
        "@type": "Answer",
        "text": "Correct across the tests you are reporting, which in an enrichment table means across gene sets. Per-gene correction happens in the differential expression step, and the pathway-level family then needs its own correction across the tested gene sets."
      }
    },
    {
      "@type": "Question",
      "name": "Is an FDR of 0.05 the same as a p-value of 0.05?",
      "acceptedAnswer": {
        "@type": "Answer",
        "text": "No. A raw p-value of 0.05 refers to a single test, while an FDR-adjusted value of 0.05 means the expected proportion of false positives among the gene sets you call significant is 5 percent. The two quantities answer different questions and adjusted values are almost always larger."
      }
    },
    {
      "@type": "Question",
      "name": "What is the best multiple testing correction for enrichment analysis?",
      "acceptedAnswer": {
        "@type": "Answer",
        "text": "For routine over-representation analysis across many gene sets, Benjamini-Hochberg FDR is the best default. For GSEA on a ranked list with adequate sample sizes, use the permutation-based FDR the method provides. For small pre-specified panels or confirmatory claims, Bonferroni or Holm family-wise control is stronger."
      }
    },
    {
      "@type": "Question",
      "name": "Why do clusterProfiler and other tools give me different adjusted p-values for the same gene list?",
      "acceptedAnswer": {
        "@type": "Answer",
        "text": "Differences usually come from default correction methods and cutoffs, the universe or background each tool uses, gene identifier mapping, and whether the tool tests over-representation or ranks. Align defaults, force a common universe and identifier type, and rerun before concluding one tool is wrong."
      }
    },
    {
      "@type": "Question",
      "name": "Can I report raw p-values alongside adjusted p-values in my enrichment table?",
      "acceptedAnswer": {
        "@type": "Answer",
        "text": "Yes. Most published tables show both, with significance claims and figures keyed off the adjusted value and raw values shown for context and labeled clearly."
      }
    },
    {
      "@type": "Question",
      "name": "My adjusted p-values come back as exactly 1.0. What does that mean?",
      "acceptedAnswer": {
        "@type": "Answer",
        "text": "Benjamini-Hochberg adjusted p-values are capped at 1, so weak tests far down the ranking are floored at 1.0. It means the test carried essentially no evidence after accounting for the number of tests, not that the calculation failed."
      }
    }
  ]
}
</script>

## Sources

1. [A pathway-level analysis reveals oxidative and inflammatory vascular stress signatures associated with VEGF inhibition.](https://pubmed.ncbi.nlm.nih.gov/42982211/)
2. [Genomic and Transcriptomic Analysis of Nine Resurrection Plant Species: Conserved Regulation Versus Structural Paralog Capacity in Desiccation Tolerance.](https://pubmed.ncbi.nlm.nih.gov/42938318/)
3. [Interpretable machine-learning prediction of PSEN1 missense variant pathogenicity based on multi-omics enrichment in six core Alzheimer's disease genes.](https://pubmed.ncbi.nlm.nih.gov/41508098/)
4. [Circulating microRNA expression profiles as biomarkers of disease activity in multiple sclerosis.](https://pubmed.ncbi.nlm.nih.gov/42499404/)
5. [Association of Lung Tissue MicroRNA Expression with Quantitative Emphysema in COPD: An Exploratory Transcriptomic Analysis.](https://pubmed.ncbi.nlm.nih.gov/42948330/)
6. [Evaluation of Two Combined Analytical Approaches to Assess the Impact of Rare Variants on Severe COVID-19.](https://pubmed.ncbi.nlm.nih.gov/42941921/)
7. [Comparative evaluation of gene selection approaches in transcriptomics: bias correction and visualization with TransPro.](https://pubmed.ncbi.nlm.nih.gov/42148814/)
8. [Analysis of Metabolic Differences Between Marselan and Cabernet Sauvignon Grape Berries from the Manas Wine-Producing Region of Xinjiang Based on Untargeted Metabolomics.](https://pubmed.ncbi.nlm.nih.gov/42883643/)
9. [Integrative multi-omics analysis of CHMP4C in colorectal cancer: associations with Genomic instability and the tumor microenvironment.](https://pubmed.ncbi.nlm.nih.gov/43007511/)
10. [Cross-scale genetic analysis prioritizes shared regulatory variation, cell-type contributions, co-expression programs, and exon-level patterns in atopic dermatitis and asthma.](https://pubmed.ncbi.nlm.nih.gov/42870119/)
11. [[Systematic evaluation of the correlation between serum metabolites and tinnitus].](https://pubmed.ncbi.nlm.nih.gov/39472117/)
12. [Rigorous validation of graph-based network analysis reappraises biologically coherent ADHD-associated transcriptomic modules in peripheral blood.](https://pubmed.ncbi.nlm.nih.gov/42832525/)