5  Enrichment foundations: ORA and GSEA

This chapter answers the first statistical decision: should a thresholded gene list be tested with over-representation analysis (ORA), or should the complete ranking be analyzed with GSEA? The examples retain the original algorithm labels so historical links remain meaningful.

5.1 Over-Representation Analysis (ORA)

Over Representation Analysis (ORA) (Boyle et al. 2004) is a widely used approach to determine whether known biological functions or processes are over-represented (= enriched) in an experimentally-derived gene list, e.g. a list of differentially expressed genes (DEGs).

enrichit implements ORA using the hypergeometric distribution (one-sided Fisher’s exact test). The p-value is calculated as the probability of observing at least \(k\) genes from the specific gene set in the selected list of \(n\) genes, given a background population (universe) of \(N\) genes containing \(M\) genes from that set.

\[ p = 1 - \sum_{i=0}^{k-1} \frac{\binom{M}{i} \binom{N-M}{n-i}}{\binom{N}{n}} \]

Example: Suppose we have 17,980 genes detected in a Microarray study and 57 genes were differentially expressed. From these, 2641 are annotated to gene set of interest. Among the differentially expressed genes, 28 belong to the gene set1.

d <- data.frame(DE_gene=c(28, 29), non_DE_gene=c(2613, 15310))
row.names(d) <- c("In_category", "not_in_category")
d

The orientation of the contingency table matters here. fisher.test() with alternative = "greater" tests whether the top-left cell is larger than expected under independence, so the top-left cell must hold the observed overlap — the 28 differentially expressed genes that belong to the gene set. Written this way, the alternative hypothesis reads exactly like our question: are the DE genes over-represented in the gene set?

Whether the gene set of interest is significantly over-represented in the differentially expressed genes can then be assessed using a hypergeometric distribution. This corresponds to a one-sided version of Fisher’s exact test.

fisher.test(d, alternative = "greater")

universe must be a character vector. It is matched against gene IDs as strings, so a numeric background (for example universe = as.numeric(ids)) matches nothing and is discarded — the analysis then silently falls back to the default background, which shifts every p-value with no sign that the argument was ignored.

Since clusterProfiler 4.21.2 you are told: the enrichment functions warn that they will ignore a non-character universe, and compareCluster() raises that warning itself — it wraps each per-cluster call in suppressMessages(), which previously swallowed the message completely. Pass the IDs as character:

universe = as.character(background_ids)

While mathematically sound, this classic model makes a naive assumption: it assumes that every gene has an equal probability of being detected. We will see how to fix this systematic bias later in the Weighted Enrichment Analysis section.

5.2 Gene Set Enrichment Analysis (GSEA)

Unlike ORA, which requires an arbitrary threshold to define “significant” genes and ignores the rest, GSEA uses the entire ranked list of genes. It calculates an Enrichment Score (ES) that reflects the degree to which a gene set is over-represented at the top or bottom of a ranked list.

There are three key elements of the GSEA method:

  • Calculation of an Enrichment Score.
    • The enrichment score (ES) represents the degree to which a set S is over-represented at the top or bottom of the ranked list L. The score is calculated by walking down the list L, increasing a running-sum statistic when we encounter a gene in S and decreasing when it is not encountered. The magnitude of the increment depends on the gene statistics (e.g., correlation of the gene with phenotype). The ES is the maximum deviation from zero encountered in the random walk; it corresponds to a weighted Kolmogorov-Smirnov(KS)-like statistic (Subramanian et al. 2005).
  • Estimation of Significance Level of ES.
    • The p-value of the ES is calculated using a permutation test. Specifically, we permute the gene labels of the gene list L and recompute the ES of the gene set for the permutated data, which generates a null distribution for the ES. The p-value of the observed ES is then calculated relative to this null distribution.
  • Adjustment for Multiple Hypothesis Testing.
    • When the entire gene sets are evaluated, the estimated significance level is adjusted to account for multiple hypothesis testing and also q-values are calculated for FDR control.

We implemented the GSEA algorithm proposed by Subramanian (Subramanian et al. 2005). enrichit offers a fast C++ implementation of GSEA. The default method (method = "multilevel") uses an adaptive multi-level splitting Monte Carlo approach (derived from fgsea (Korotkevich et al. 2019)) to estimate low p-values efficiently with high accuracy.

library(enrichit)

# Generate synthetic ranked gene list
set.seed(42)
geneList <- sort(rnorm(1000), decreasing = TRUE)
names(geneList) <- paste0("Gene", 1:1000)

gene_sets <- list(
  PathwayTop = names(geneList)[1:50],
  PathwayBottom = names(geneList)[951:1000]
)

# Run GSEA using the multilevel method
gsea_result <- gsea(
  geneList = geneList,
  gene_sets = gene_sets,
  method = "multilevel"
)

5.2.1 Leading edge analysis and core enriched genes

Leading edge analysis reports Tags to indicate the percentage of genes contributing to the enrichment score, List to indicate where in the list the enrichment score is attained and Signal for enrichment signal strength.

It would also be very interesting to get the core enriched genes that contribute to the enrichment. Our packages (clusterProfiler, DOSE, meshes and ReactomePA) support leading edge analysis and report core enriched genes in GSEA analysis.

5.2.2 One-sided or two-sided: the scoreType argument

scoreType controls which end of the ranked list the enrichment score is allowed to look at. It takes three values, and the choice is not cosmetic — it decides which gene sets are even visible to the test:

scoreType ES is Detects
"std" (default) the maximum deviation from zero in either direction gene sets enriched at the top and at the bottom
"pos" the maximum deviation above zero only gene sets enriched at the top only
"neg" the maximum deviation below zero only gene sets enriched at the bottom only

With the usual signed metric (a log fold change or a signed correlation), the default "std" is what you want: it reports both up- and down-regulated sets. Running the three settings on the same mixed-sign ranking makes the difference obvious — "pos" misses the bottom-enriched set entirely (its ES is 0 and its p-value is 1), and "neg" misses the top-enriched one.

## gene sets enriched at the top and at the bottom of the same ranking
gsea(geneList, gene_sets, scoreType = "std")   # finds both
gsea(geneList, gene_sets, scoreType = "pos")   # only the top set
gsea(geneList, gene_sets, scoreType = "neg")   # only the bottom set

The case that trips people up is a metric that has no sign at all — a score where higher simply means “more strongly associated with the phenotype”, such as an evidence score, a correlation magnitude, or a rank-based statistic that runs from 1 to \(n\). Here every value is positive, so the “bottom” of the list is not a meaningful direction and "std" is the wrong model. Use scoreType = "pos":

## all-positive evidence scores: higher = more associated
gsea(evidence_score, gene_sets, scoreType = "pos")

The package warns you when it can detect this situation (All values in the stats vector are greater than zero and scoreType is "std", maybe you should switch to scoreType = "pos"), but the warning is easy to ignore — and it only fires when every value is positive, so it will not catch a metric that is positive except for a handful of zeros. Note also that the normalized enrichment score (NES) is computed against a null that depends on scoreType, so comparing NES values obtained with different settings is not meaningful.

5.2.3 Choosing between ORA and GSEA

A question that comes up constantly in teaching: given my experiment, which one should I run? The deciding factor is mostly how many genes pass your significance threshold — not personal preference, and not the analysis software you happen to have installed.

Situation Recommended Why
Many DEGs (rule of thumb: dozens and up) ORA The gene list is long enough for the hypergeometric test to have real power, and the result is easy to read term by term.
Few DEGs (say, under ~20–30) GSEA ORA has almost no power on a short list: with \(n\) small, only a very strong overlap can reach significance. GSEA ranks the whole list, so a modest but coordinated shift across a gene set is still detectable.
No DEG passes the cutoff at all GSEA ORA has literally nothing to test. GSEA only needs a ranking, so it still works when the threshold yields an empty list.
You do not want to commit to an arbitrary cutoff GSEA No threshold is involved anywhere in the method.
You want a per-term statement of the form “this set is over-represented in my DE genes” ORA That is exactly the question the hypergeometric test answers.

The failure mode to avoid is running ORA on a handful of genes. With, say, ten DEGs the test can only detect very large overlaps, so it typically returns nothing — and it is easy to misread that null result as “no biology here” rather than “not enough data for this test”. In that regime GSEA (or the network-based and weighted methods introduced below) is the appropriate tool.

The two methods are not mutually exclusive, and reporting both is common. They answer slightly different questions: ORA asks whether your thresholded genes hit a set more often than chance, while GSEA asks whether a set is shifted toward one end of the full ranking.

5.3 Next steps

References

Boyle, Elizabeth I, Shuai Weng, Jeremy Gollub, et al. 2004. “GO::TermFinder–open Source Software for Accessing Gene Ontology Information and Finding Significantly Enriched Gene Ontology Terms Associated with a List of Genes.” Bioinformatics (Oxford, England) 20 (18): 3710–15. https://doi.org/10.1093/bioinformatics/bth456.
Korotkevich, Gennady, Vladimir Sukhov, and Alexey Sergushichev. 2019. “Fast Gene Set Enrichment Analysis.” bioRxiv, October 22, 060012. https://doi.org/10.1101/060012.
Subramanian, Aravind, Pablo Tamayo, Vamsi K. Mootha, et al. 2005. “Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-Wide Expression Profiles.” Proceedings of the National Academy of Sciences of the United States of America 102 (43): 15545–50. https://doi.org/10.1073/pnas.0506580102.

  1. example adopted from https://guangchuangyu.github.io/2012/04/enrichment-analysis/↩︎