Functional Enrichment

Turn a DE gene list into pathways and GO terms with clusterProfiler.

Overview

Problem. Translate a long DE list into interpretable pathways.

Learning goals

• GSEA (full ranking) vs ORA (thresholded)
• Enrichment generates hypotheses, not conclusions

Figures

Tutorial

Translate differential expression results into biological insights using GSEA and ORA.

When to Use This Skill

Use this skill after completing differential expression analysis to identify enriched pathways and biological processes.

Two complementary methods: GSEA (primary, uses all ranked genes, detects coordinated changes) and ORA (secondary, uses significant gene list, validates GSEA). Default recommendation: Run GSEA unless user specifically requests ORA or has only a gene list (no fold changes).

See references/gsea_ora_comparison.md for detailed method comparison.

Quick Start (Example Data)

Test this skill with real DE results in ~2 minutes:

# Load example DE results from airway dataset (dexamethasone treatment)
source("scripts/load_example_data.R")
de_results <- load_airway_de_results()  # Auto-installs packages (~1-2 min, ~40MB)

# Load required packages and scripts
library(clusterProfiler)
library(msigdbr)
source("scripts/prepare_gene_lists.R")
source("scripts/get_msigdb_genesets.R")
source("scripts/run_gsea.R")
source("scripts/generate_plots.R")
source("scripts/export_results.R")

# Run GSEA workflow
ranked_genes <- create_ranked_list(de_results)
term2gene <- get_msigdb_genesets("human", c("H"))  # Hallmark pathways only for speed
gsea_result <- run_gsea(ranked_genes, term2gene, n_perm = 1000)
generate_all_plots(gsea_result)
export_all(gsea_result, ranked_genes, output_prefix = "quick_test")

For your own data: Replace example data loading with your DE results file (see Inputs section).

Installation

# Install Bioconductor packages
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install(c("clusterProfiler", "enrichplot", "org.Hs.eg.db"))

# Install CRAN packages
install.packages(c("msigdbr", "ggplot2", "ggrepel", "dplyr"))

Software Requirements

Database Licenses

⚠️ KEGG requires commercial license for commercial use. MSigDB and GO use CC-BY 4.0 (commercial use permitted with attribution). For commercial applications, use Hallmark, Reactome, or GO databases, or obtain appropriate KEGG license.

Inputs

Outputs

Clarification Questions

⚠️ CRITICAL: Always ask question #1 first to check if user has provided input files before proceeding with analysis.

2. Method: GSEA (default with full DE results), ORA (gene list only), or Both (validation)
3. Databases: Hallmark + KEGG (default exploratory), GO:BP (detailed processes), Reactome (alternative to KEGG), C6/C7 (cancer/immunology)
4. Species: Human (default) or Mouse

See references/decision-guide.md for detailed decision guidance.

Standard Workflow

🚨 MANDATORY: USE SCRIPTS EXACTLY AS SHOWN - DO NOT WRITE INLINE CODE 🚨

CRITICAL: Use relative paths (scripts/). DO NOT construct absolute paths.

These scripts provide robust, validated enrichment analysis. Use them as-is for standard GSEA workflows.

Step 1 - Load data:

# Load all required scripts
source("scripts/load_de_results.R")
source("scripts/prepare_gene_lists.R")
source("scripts/get_msigdb_genesets.R")
source("scripts/run_gsea.R")
source("scripts/generate_plots.R")
source("scripts/export_results.R")

# Load your DE results
de_results <- load_de_results("your_de_results.csv")  # Replace with your file path

# Prepare inputs
ranked_genes <- create_ranked_list(de_results)
term2gene <- get_msigdb_genesets("human", c("H", "C2:CP:KEGG"))

Step 2 - Run analysis:

gsea_result <- run_gsea(ranked_genes, term2gene)

Step 3 - Generate visualizations:

generate_all_plots(gsea_result)

The script handles PNG + SVG export with graceful fallback for SVG dependencies.

Step 4 - Export results:

export_all(gsea_result, ranked_genes, output_prefix = "enrichment")

❌ IF YOU DON'T SEE THESE: You wrote inline code. Stop and use the scripts.

NEVER skip directly to writing inline code without trying the script first.

Decision Guide

Key decisions during analysis:

1. Method selection (before analysis): GSEA only (default, full DE results), ORA only (gene list without fold changes), or Both (validation)

2. Database selection (Step 3): Hallmark + KEGG (default exploratory), GO:BP (detailed processes), Reactome (KEGG alternative, commercial-friendly), or C6/C7 (cancer/immunology)

3. Ranking metric (Step 2): Test statistic (best, balances effect + significance), signed -log10(pvalue) (good alternative), or log2FC only (not recommended)

For detailed decision guidance with options and recommendations, see references/decision-guide.md.

Common Patterns

For all patterns and advanced examples, see references/method-reference.md.

Pattern 1: Quick GSEA with Default Settings

# Load required scripts
source("scripts/load_de_results.R")
source("scripts/prepare_gene_lists.R")
source("scripts/get_msigdb_genesets.R")
source("scripts/run_gsea.R")
source("scripts/generate_plots.R")
source("scripts/export_results.R")

# Load, prepare, and run GSEA
de_results <- load_de_results("de_results.csv")
ranked_genes <- create_ranked_list(de_results)
term2gene <- get_msigdb_genesets("human", c("H", "C2:CP:KEGG"))
gsea_result <- run_gsea(ranked_genes, term2gene)
generate_all_plots(gsea_result)
export_all(gsea_result, ranked_genes, output_prefix = "enrichment")

Pattern 2: GSEA + ORA for Validation

# Load required scripts
source("scripts/load_de_results.R")
source("scripts/prepare_gene_lists.R")
source("scripts/get_msigdb_genesets.R")
source("scripts/run_gsea.R")
source("scripts/run_ora.R")
source("scripts/generate_plots.R")
source("scripts/export_results.R")

# Prepare inputs and run both methods
de_results <- load_de_results("de_results.csv")
ranked_genes <- create_ranked_list(de_results)
sig_genes <- filter_significant_genes(de_results)
term2gene <- get_msigdb_genesets("human", c("H", "C2:CP:KEGG"))
gsea_result <- run_gsea(ranked_genes, term2gene)
ora_up <- run_ora(sig_genes$up, term2gene, sig_genes$background, direction = "upregulated")
ora_down <- run_ora(sig_genes$down, term2gene, sig_genes$background, direction = "downregulated")
generate_all_plots(gsea_result)
generate_ora_barplot(ora_up, "Upregulated")
generate_ora_barplot(ora_down, "Downregulated")
export_all(gsea_result, ora_up, ora_down, ranked_genes, output_prefix = "enrichment")

For results interpretation guidance, see references/interpretation_guidelines.md.

Common Errors

Suggested Next Steps

After completing functional enrichment analysis:

1. Interpret and Validate Results

Check if enriched pathways align with experimental hypothesis and are biologically coherent. Examine leading edge genes and compare with literature. See references/interpretation_guidelines.md for detailed guidance.

2. Create Additional Visualizations

Generate optional advanced plots: enrichment maps (pathway networks), concept networks (gene-pathway connections), or ridge plots (expression distributions). Use scripts/generate_plots.R.

3. Deep Dive on Key Pathways

Examine leading edge genes in top pathways, check individual gene expression in original DE results, and look up pathway diagrams in KEGG/Reactome databases.

4. Export for Publications

Use CSV results for supplementary tables, SVG plots for figures (300+ DPI), and include enrichment summary in results section. Document database versions, parameters, and software versions in methods.

Related Skills

References

Primary Citations

1. clusterProfiler: Yu G, et al. (2012) clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. DOI: 10.1089/omi.2011.0118

2. GSEA method: Subramanian A, et al. (2005) Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. PNAS. DOI: 10.1073/pnas.0506580102

3. MSigDB: Liberzon A, et al. (2015) The Molecular Signatures Database Hallmark Gene Set Collection. Cell Systems. DOI: 10.1016/j.cels.2015.12.004

4. Method benchmarking: Geistlinger L, et al. (2021) Toward a gold standard for benchmarking gene set enrichment analysis. Brief Bioinform. DOI: 10.1093/bib/bbz158

Additional Resources

Code preview

scripts/export_results.R

#' Export Enrichment Results
#'
#' Functions to export GSEA and ORA results to CSV files and generate summary reports.


#' Export Results to CSV
#'
#' Saves GSEA and ORA results to CSV files.
#'
#' @param gsea_result enrichResult object from run_gsea()
#' @param ora_up enrichResult object from run_ora() for upregulated genes
#' @param ora_down enrichResult object from run_ora() for downregulated genes
#' @param output_prefix Prefix for output filenames (default: "enrichment")
#' @export
#'
#' @examples
#' export_results(gsea_result, ora_up, ora_down, output_prefix = "my_analysis")

export_results <- function(gsea_result = NULL, ora_up = NULL, ora_down = NULL, output_prefix = "enrichment") {

  message("\n=== Exporting Results ===")

  # GSEA results
  if (!is.null(gsea_result) && nrow(gsea_result@result) > 0) {
    gsea_df <- as.data.frame(gsea_result@result)
    gsea_file <- paste0(output_prefix, "_gsea_results.csv")
    write.csv(gsea_df, gsea_file, row.names = FALSE)
    message(sprintf("  Saved: %s (%d rows)", gsea_file, nrow(gsea_df)))
  }

  # ORA results (up)
  if (!is.null(ora_up) && nrow(ora_up@result) > 0) {
    ora_up_df <- as.data.frame(ora_up@result)
    ora_up_df$direction <- "upregulated"
    ora_up_file <- paste0(output_prefix, "_ora_up_results.csv")
    write.csv(ora_up_df, ora_up_file, row.names = FALSE)
    message(sprintf("  Saved: %s (%d rows)", ora_up_file, nrow(ora_up_df)))
  }

  # ORA results (down)
  if (!is.null(ora_down) && nrow(ora_down@result) > 0) {
    ora_down_df <- as.data.frame(ora_down@result)
    ora_down_df$direction <- "downregulated"
    ora_down_file <- paste0(output_prefix, "_ora_down_results.csv")
    write.csv(ora_down_df, ora_down_file, row.names = FALSE)
    message(sprintf("  Saved: %s (%d rows)", ora_down_file, nrow(ora_down_df)))
  }
}


#' Generate Summary Report
#'
#' Creates a markdown summary of top enrichment findings.
#'
#' @param gsea_result enrichResult object from run_gsea()
#' @param ora_up enrichResult object from run_ora() for upregulated genes
#' @param ora_down enrichResult object from run_ora() for downregulated genes
#' @param output_file Output filename (default: "enrichment_summary.md")
#' @export
#'
#' @examples
#' generate_summary(gsea_result, ora_up, ora_down, output_file = "my_summary.md")

#' Export All Results (export_all wrapper)
#'
#' Comprehensive export function that saves all enrichment results including
#' CSV files, markdown summary, and RDS objects for downstream analysis.
#'
#' This is the recommended function to call for complete result export.
#'
#' @param gsea_result enrichResult object from run_gsea()
#' @param ora_up enrichResult object from run_ora() for upregulated genes (optional)
#' @param ora_down enrichResult object from run_ora() for downregulated genes (optional)
#' @param ranked_genes Named numeric vector of ranked genes (optional, for downstream use)
#' @param output_prefix Prefix for output filenames (default: "enrichment")
#' @export
#'
#' @examples
#' export_all(gsea_result, output_prefix = "my_analysis")
#' export_all(gsea_result, ora_up, ora_down, ranked_genes, output_prefix = "my_analysis")

scripts/generate_plots.R

#' Visualization Functions for Enrichment Results
#'
#' This script contains functions to generate standard plots for GSEA and ORA results.

library(enrichplot)
library(ggplot2)
library(dplyr)

# Try to load svglite for high-quality SVG (optional)
.has_svglite <- requireNamespace("svglite", quietly = TRUE)
if (.has_svglite) {
  library(svglite)
}

#' Save Plot in Both PNG and SVG Formats
#'
#' Internal helper to save plots with graceful fallback for SVG.
#'
#' @param plot ggplot object
#' @param base_path Base filename (without extension)
#' @param width Width in inches
#' @param height Height in inches
#' @param dpi DPI for PNG output
.save_plot <- function(plot, base_path, width = 8, height = 6, dpi = 300) {
  # Always save PNG
  png_path <- sub("\\.(svg|png)$", ".png", base_path)
  ggsave(png_path, plot = plot, width = width, height = height, dpi = dpi, device = "png")
  cat("   Saved:", png_path, "\n")

  # Always try SVG - try ggsave first, fall back to svg() device
  svg_path <- sub("\\.(svg|png)$", ".svg", base_path)
  tryCatch({
    ggsave(svg_path, plot = plot, width = width, height = height, device = "svg")
    cat("   Saved:", svg_path, "\n")
  }, error = function(e) {
    # If ggsave fails, try base R svg() device directly
    tryCatch({
      svg(svg_path, width = width, height = height)
      print(plot)
      dev.off()
      cat("   Saved:", svg_path, "\n")
    }, error = function(e2) {
      cat("   (SVG export failed)\n")
    })
  })
}

#' Generate GSEA Dot Plot (PRIMARY - always generate)
#'
#' Creates a dot plot showing activated vs suppressed pathways from GSEA.
#'
#' @param gsea_result enrichResult object from run_gsea()
#' @param top_n Number of top pathways to show (default: 20)
#' @param output_file Output filename (default: "gsea_dotplot.svg")
#' @return ggplot object
#' @export

generate_gsea_dotplot <- function(gsea_result, top_n = 20, output_file = "gsea_dotplot.svg") {
  if (is.null(gsea_result) || nrow(gsea_result@result) == 0) {
    message("No GSEA results to plot")
    return(NULL)
  }

  p <- dotplot(gsea_result, showCategory = top_n, split = ".sign") +
    facet_grid(~.sign) +
    theme_bw() +
    theme(
      axis.text.y = element_text(size = 8),
      strip.text = element_text(size = 10, face = "bold")
    ) +
    ggtitle("GSEA Results: Activated vs Suppressed Pathways")

  .save_plot(p, output_file, width = 12, height = 8)

  return(p)
}


#' Generate GSEA Running Score Plot
#'

scripts/get_msigdb_genesets.R

#' Get MSigDB Gene Sets
#'
#' Retrieves gene sets from MSigDB for use in enrichment analysis.
#' Supports multiple databases and species.
#'
#' @param species Species name: "human" or "mouse" (default: "human")
#' @param categories Vector of category names to include (default: c("H"))
#'   Options: "H" (Hallmark), "KEGG", "REACTOME", "GO:BP", "GO:MF", "GO:CC"
#' @return Data frame in TERM2GENE format (columns: term, gene)
#' @export
#'
#' @examples
#' # Hallmark gene sets only (default)
#' term2gene <- get_msigdb_genesets(species = "human", categories = c("H"))
#'
#' # Hallmark + KEGG
#' term2gene <- get_msigdb_genesets(species = "human", categories = c("H", "KEGG"))
#'
#' # GO Biological Process
#' term2gene <- get_msigdb_genesets(species = "human", categories = c("GO:BP"))

library(msigdbr)

get_msigdb_genesets <- function(species = "human", categories = c("H")) {
  # Map species
  species_name <- ifelse(species == "human", "Homo sapiens", "Mus musculus")

  # Category mapping (updated for msigdbr >= 10.0.0)
  cat_map <- list(
    "H" = list(collection = "H", subcollection = NULL),           # Hallmark
    "KEGG" = list(collection = "C2", subcollection = "CP:KEGG"),  # KEGG
    "REACTOME" = list(collection = "C2", subcollection = "CP:REACTOME"),
    "GO:BP" = list(collection = "C5", subcollection = "GO:BP"),
    "GO:MF" = list(collection = "C5", subcollection = "GO:MF"),
    "GO:CC" = list(collection = "C5", subcollection = "GO:CC"),
    # Support for old-style subcategory format (e.g., "C2:CP:KEGG")
    "C2:CP:KEGG" = list(collection = "C2", subcollection = "CP:KEGG"),
    "C2:CP:REACTOME" = list(collection = "C2", subcollection = "CP:REACTOME"),
    "C5:GO:BP" = list(collection = "C5", subcollection = "GO:BP"),
    "C5:GO:MF" = list(collection = "C5", subcollection = "GO:MF"),
    "C5:GO:CC" = list(collection = "C5", subcollection = "GO:CC")
  )

  all_genesets <- list()

  for (cat in categories) {
    if (!cat %in% names(cat_map)) {
      warning(sprintf("Unknown category: %s", cat))
      next
    }

    params <- cat_map[[cat]]

    # Use new msigdbr API (collection/subcollection instead of category/subcategory)
    if (is.null(params$subcollection)) {
      msig <- msigdbr(species = species_name, collection = params$collection)
    } else {
      msig <- msigdbr(species = species_name, collection = params$collection,
                      subcollection = params$subcollection)
    }

    # Convert to list format for clusterProfiler
    gs <- split(msig$gene_symbol, msig$gs_name)
    all_genesets <- c(all_genesets, gs)

    message(sprintf("Loaded %d gene sets from %s", length(unique(msig$gs_name)), cat))
  }

  # Convert to TERM2GENE format
  term2gene <- do.call(rbind, lapply(names(all_genesets), function(term) {
    data.frame(term = term, gene = all_genesets[[term]], stringsAsFactors = FALSE)
  }))

  return(term2gene)
}

Companion files