Chapter 11: Gene Set Enrichment Analysis (GSEA) & clusterProfiler Workflows

1. Why Gene Set Enrichment Analysis (GSEA) Outperforms ORA

Over-Representation Analysis (ORA) relies on an arbitrary significance threshold (e.g., \(p < 0.05\) and \(|\log_2\text{FC}| \ge 1.0\)). This introduces two major weaknesses:

1. Genes just below the threshold (e.g., \(p = 0.051\)) are completely discarded.
2. It treats all selected genes equally, ignoring their actual effect sizes.

Developed by Subramanian et al. (Broad Institute), Gene Set Enrichment Analysis (GSEA) does not filter by arbitrary cutoffs. Instead, it utilizes the entire list of 20,000 genes, ranked by a signed metric combining fold change and \(p\)-value:

$$S_i = \text{sign}(\log_2\text{FC}_i) \times \left(-\log_{10} p_i\right)$$

The GSEA Running Sum Algorithm

GSEA steps through the ranked gene list from top to bottom. When a gene belongs to the pathway, the running sum score increases; when it does not, the score decreases. The Enrichment Score (ES) is the maximum deviation of this running sum from zero. Permutation testing determines statistical significance, normalized by gene set size into a Normalized Enrichment Score (NES).

2. Production GSEA Script with MSigDB Hallmark Signatures

library(clusterProfiler)
library(msigdbr) # Molecular Signatures Database access
library(enrichplot)

# 1. Prepare Full Ranked Gene Vector
# Ranking metric = sign(log2FC) * -log10(pvalue)
res_clean <- res_df[!is.na(res_df$padj) & !is.na(res_df$log2FoldChange), ]
res_clean$rank_metric <- sign(res_clean$log2FoldChange) * -log10(res_clean$pvalue + 1e-300)

# Sort descending
ranked_genes <- res_clean$rank_metric
names(ranked_genes) <- res_clean$Gene_Symbol
ranked_genes <- sort(ranked_genes, decreasing = TRUE)

# 2. Retrieve MSigDB Hallmark Gene Sets (50 Well-Defined Biological States)
h_gene_sets <- msigdbr(species = "Homo sapiens", category = "H") %>%
    dplyr::select(gs_name, gene_symbol)

# 3. Execute GSEA
gsea_results <- GSEA(
    geneList = ranked_genes,
    TERM2GENE = h_gene_sets,
    pvalueCutoff = 0.05,
    pAdjustMethod = "BH",
    verbose = FALSE
)

# 4. GSEA Visualizations
# Multi-pathway dotplot showing positive and negative NES
dotplot(gsea_results, showCategory = 10, title = "Hallmark GSEA: Enriched Pathways")

# Running score mountain plot for top pathway
gseaplot2(gsea_results, geneSetID = 1, title = gsea_results$Description[1])

# Ridgeplot showing expression distribution shifts across enriched pathways
ridgeplot(gsea_results, showCategory = 8)