GSEA Teaching Manual: HALLMARK_G2M_CHECKPOINT as an Example

GSEA and ORA: names, method and differences

GSEA stands for Gene Set Enrichment Analysis. ORA stands for Over-Representation Analysis; what people usually call "GO / KEGG enrichment analysis" is most often ORA.

ORA (for example a hypergeometric test on GO or KEGG terms) first cuts a DEG list with thresholds such as padj < 0.05 and |log2FC| > 1, then asks whether pathway members are over-represented in that list. GSEA uses no threshold. It ranks all genes by a differential metric and asks whether the gene set as a whole drifts toward one end of the ranked list.

How ORA is computed

ORA asks: are genes from a given pathway more common in the significant-gene list than random sampling would predict? It uses a hypergeometric test (equivalent to a one-sided Fisher's exact test) on this 2×2 table.

N is the number of background genes (the universe), K the number of genes in the pathway, n the number of DEGs, and k the overlap. The p-value is the probability of drawing at least k pathway members when n genes are drawn at random from N:

Typical reports show GeneRatio = k/n and BgRatio = K/N. A GeneRatio clearly above the BgRatio means the pathway is over-represented among the DEGs.

Side-by-side comparison

The Hallmark collection (MSigDB Collection H) distils thousands of overlapping MSigDB gene sets into 50 well-defined, minimally overlapping sets, each representing a specific biological state or process. Running these 50 first keeps the multiple-testing burden small and the results focused.

The standard GSEA workflow

1. Prepare expression data and groupsAn RNA-seq count matrix or a microarray expression matrix, plus two phenotype labels (for example Tumor vs Normal). The GSEA desktop app uses .gct and .cls files.
2. Compute a ranking metric r_i to build the ranked list LCommon choices are Signal-to-Noise (the desktop default), a t statistic, the DESeq2 Wald stat, or sign(log2FC) × −log10(p). log2FC alone is not recommended because it ignores variance. Sort the list by r from high to low; collapse duplicate gene symbols first.
3. Load gene sets (.gmt)For example the human Hallmark file h.all.v2024.1.Hs.symbols.gmt (the version number depends on when you download it). Usually only gene sets with 15–500 members present in the data are kept.
4. Walk the list to get the running ES and the ESMembers (hits) add a score weighted by |r_i|^p; non-members (misses) subtract an equal share. The standard weighting is p = 1 (weighted).
5. Permutation test to build a null distributionRepeat 1,000 times. With ≥ 7 samples per group use phenotype permutation (shuffle sample labels, which keeps gene–gene correlation); with fewer samples, or with Preranked, use gene_set permutation.
6. Normalize and correct: NES, nominal p, FDR qNES divides by the mean of same-sign null ES values to remove the effect of gene set size; FDR evaluates all gene sets together.
7. Leading-edge analysisFind the core subset of genes that pushes the ES to its peak.
8. Interpret and validateFilter by NES / FDR, check the shape of the plot, then pick leading-edge genes for qPCR, Western blots or functional experiments.

Example gene set: HALLMARK_G2M_CHECKPOINT

This set has 200 genes. It represents the checkpoint and division machinery between the end of DNA replication and the start of mitosis: the cyclin B–CDK1 complex, PLK1 and AURKA that drive entry into mitosis, CHEK1 that brakes the cycle on DNA damage, the CDC25 phosphatases that activate CDK1 by dephosphorylation, and the spindle assembly checkpoint genes BUB1, MAD2L1 and CDC20.

Interpretation

Statistical thresholds

Four common curve shapes

Interpretation checklist

• Direction: the sign of NES is relative to how the list is ranked. Check whether r is “treatment − control” or the reverse; otherwise up and down flip completely.
• Strength and confidence: discuss only |NES| > 1.5 with FDR < 0.25; for main conclusions in a paper, aim for FDR < 0.05.
• Shape: is the peak sharp and near the edge, and is the barcode really packed at one end? If the curve is flat and the barcode evenly spread, be cautious even if the p value looks good.
• Leading-edge size: if only 10 of 200 members are in the leading edge, a few genes drive the signal; check whether those genes are outliers themselves.
• Redundancy: when G2M_CHECKPOINT, E2F_TARGETS and MITOTIC_SPINDLE are all significant, describe them as one “proliferation” story.
• Sample size and permutation type: with < 7 samples per group only gene_set permutation is possible; it ignores gene–gene correlation, so the FDR is optimistic and conclusions should be worded cautiously.

In practice: R (fgsea + msigdbr)

# DESeq2 results in res (Tumor vs Normal), with a gene symbol column `symbol` added
library(msigdbr); library(fgsea); library(dplyr)

# From msigdbr 10 the full data live in a separate package, msigdbdf; install it once
h  <- msigdbr(species = "Homo sapiens", collection = "H")   # older versions: category = "H"
pw <- split(h$gene_symbol, h$gs_name)

rk <- res |> as.data.frame() |> filter(!is.na(stat), !is.na(symbol)) |>
      group_by(symbol) |> summarise(stat = mean(stat))      # collapse duplicate genes
ranks <- sort(setNames(rk$stat, rk$symbol), decreasing = TRUE)

set.seed(42)
fg <- fgsea(pw, ranks, minSize = 15, maxSize = 500)
fg |> arrange(padj) |> select(pathway, NES, padj, size) |> head(10)

plotEnrichment(pw[["HALLMARK_G2M_CHECKPOINT"]], ranks) +
  ggplot2::labs(title = "HALLMARK_G2M_CHECKPOINT")
fg$leadingEdge[[which(fg$pathway == "HALLMARK_G2M_CHECKPOINT")]]

fgsea’s padj is a Benjamini–Hochberg adjusted p value. It is computed differently from the FDR q of the GSEA desktop app but serves the same purpose. In the desktop app, suggested settings are: Permutations 1000, Enrichment statistic weighted, Metric Signal2Noise, Max/Min size 500/15.

Common misconceptions