GPTomics/bio-pathway-gsea

Version Compatibility

Reference examples tested with: clusterProfiler 4.18.4+, org.Hs.eg.db 3.22+, msigdbr 26+, fgsea 1.36+.

Before using code patterns, verify installed versions match. If versions differ:

• R: packageVersion('<pkg>') then ?function_name to verify parameters

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

gseKEGG queries the live KEGG REST API, so the same code returns different results as KEGG updates; pin the run date. gseGO/gsePathway and MSigDB GSEA use local annotation (org.*.eg.db, reactome.db, msigdbr) and are reproducible given the package version. The single source of truth for versions is this block, not headings.

Gene Set Enrichment Analysis (GSEA)

"Which pathways shift coordinately across my full ranked gene list, with no cutoff?" -> Walk a weighted running-sum down the genome-wide ranking and test whether each gene set piles up at one END - because that score reports the structure of YOUR ranking, so the ranking metric and the permutation type, not the gene sets, decide the result.

• R: gseGO(geneList, OrgDb, ont), gseKEGG(geneList, organism), GSEA(geneList, TERM2GENE)

Scope: threshold-free Functional Class Scoring (FCS) of a RANKED vector - the running-sum ES, the ranking-metric choice, the permutation null, NES/FDR, the leading edge, and per-sample ssGSEA/GSVA scores. A pre-selected unranked gene LIST -> go-enrichment (ORA). The ranking statistic source -> differential-expression/de-results. KEGG/Reactome/WikiPathways database semantics -> kegg-pathways, reactome-pathways, wikipathways. Plots -> enrichment-visualization.

The Single Most Important Modern Insight -- A GSEA Result Is a Deterministic Function of Three Usually-Implicit Choices: the Ranking Statistic, the Weight Exponent p, and the Permuted Labels

GSEA is not a discovery about biology - it is the running-sum's report on the ranking it was handed. The weighted enrichment score (Subramanian 2005; the default exponent=1 weights each hit by the gene's statistic magnitude) asks exactly one question: do the members of a set pile up at one END of YOUR ranking. So the result is fixed by three decisions tutorials usually leave silent, and Wijesooriya 2022 found most published GSEA papers report none of them.

1. The ranking statistic IS the experiment. Rank by a signed, variance-calibrated metric (DESeq2 stat, limma moderated t) - sign gives direction, variance-calibration sinks noisy low-information genes to the middle. Ranking by a RAW p-value erases the sign, so up- and down-regulated genes collapse together and NES becomes uninterpretable. Ranking by bare log2FC lets a handful of low-count genes with huge unstable fold changes hijack the leading edge. A bad ranking is faithfully reported as a ranking artifact.
2. The permutation type sets validity. Phenotype (sample-label) permutation preserves the gene-gene correlation that co-regulated pathways are made of - it is the gold standard for type-I error control because it preserves that correlation, but needs the expression matrix and adequate n per group (~>=7); below that n its validity degrades. clusterProfiler/fgsea preranked are FORCED into GENE permutation, which treats genes as independent, destroys that correlation, and is ANTI-CONSERVATIVE: it returns "significant" pathways that are nothing but co-expression. Accept it, report it, and prefer CAMERA when the design matrix is available (full competitive/self-contained theory in the category README).
3. The honesty is bounded by reporting. Log the ranking metric, the exponent p, the permutation type, the gene-set collection and its version/date, the size filters, and the multiple-testing method. Without those, the result is unfalsifiable.

Tool Taxonomy

| Method | Citation | Mechanism / role | When | |--------|----------|------------------|------| | Preranked GSEA (gseGO/gseKEGG/GSEA) | Subramanian 2005 PNAS 102:15545; Mootha 2003 Nat Genet 34:267 | weighted running-sum ES over a ranked vector; gene-permutation null | the common case: a ranked statistic for all genes, no matrix | | fgsea engine | Korotkevich 2021 bioRxiv 060012 (preprint) | fgseaMultilevel; resolves tiny p accurately down to eps | the engine under by='fgsea' (default); what gives sub-1/nperm p-values | | Phenotype-permutation GSEA | Subramanian 2005 PNAS 102:15545 | shuffles sample labels; preserves gene-gene correlation | matrix + phenotype + ~>=7/group; the gold-standard competitive test | | CAMERA | Wu & Smyth 2012 NAR 40:e133 | competitive, VIF-corrects inter-gene correlation analytically | matrix + design; want a correlation-honest competitive test | | ROAST / fry | Wu 2010 Bioinformatics 26:2176 | self-contained rotation test; valid at any n | matrix + design, tiny n, "is the set DE at all" | | ssGSEA | Barbie 2009 Nature 462:108 | per-sample rank-based enrichment score | a per-sample pathway-activity matrix (no contrast test) | | GSVA | Hanzelmann 2013 BMC Bioinformatics 14:7 | unsupervised per-sample, per-set kernel/CDF score | per-sample features for clustering/survival/ML |

Decision Tree by Scenario

| Scenario | Recommended | Why | |----------|-------------|-----| | Ranked statistic for ALL genes, cutoff would be arbitrary | preranked gseGO/gseKEGG/GSEA | threshold-free; the common case (report gene-permutation) | | Pre-selected unranked list (module, GWAS hits, screen) | ORA -> go-enrichment | no genome-wide ranking exists | | Function annotation, broad GO coverage | gseGO | GO is the broadest local resource | | Metabolic / signaling pathways | gseKEGG -> kegg-pathways | KEGG maps (live DB) | | Reaction-level, reproducible offline | gsePathway -> reactome-pathways | local reactome.db | | Curated MSigDB hallmark / C2 / C5 | GSEA(TERM2GENE) + msigdbr | generic-input GSEA on any collection | | Matrix + design, want competitive + correlation-honest | limma::camera | VIF-corrects the inter-gene correlation gene-permutation ignores | | Matrix + design, tiny n / covariates, "is set DE at all" | limma::roast/fry | self-contained rotation, valid at any n | | Per-sample pathway-activity matrix for clustering/ML | ssGSEA / GSVA | scores each sample, not a contrast test | | The DE statistic / ranking itself | -> differential-expression/de-results | upstream, not enrichment |

Build the Ranked Vector (the Caller-Owned Step)

Goal: Turn a DE table into the named numeric vector sorted strictly decreasing that every preranked function requires, ranked by a signed variance-calibrated metric.

Approach: Pick the ranking metric to match the DE tool, name the vector by gene ID, drop NAs, deduplicate to one statistic per gene, and sort decreasing. The DE mechanics and the $padj/$adj.P.Val column conventions live in differential-expression/de-results.

| DE source | Ranking metric | Column | Why | |-----------|----------------|--------|-----| | DESeq2 | Wald statistic | stat | signed + variance-calibrated; best single choice for RNA-seq | | limma / voom | moderated t-statistic | t | empirical-Bayes shrinkage borrows variance; signed | | edgeR (QL) | sign(logFC) * -log10(PValue) | derived | no single signed statistic column | | any tool, last resort | log2 fold change | log2FoldChange | magnitude only; noisy for low-count genes |

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

de <- read.csv('de_results.csv')          # DE list source: differential-expression/de-results
gene_list <- de$stat                       # DESeq2 Wald stat: signed + variance-calibrated
names(gene_list) <- de$entrez_id
gene_list <- gene_list[!is.na(gene_list)]
gene_list <- gene_list[!duplicated(names(gene_list))]   # one statistic per gene; duplicates double-count hits
gene_list <- sort(gene_list, decreasing = TRUE)         # REQUIRED: unsorted input silently mis-ranks

Ranking by sign(log2FC) * -log10(pmax(pvalue, 1e-300)) (for edgeR, or when a Wald stat is unavailable) preserves direction and clamps p==0 from going to Inf. Never rank by raw p-value alone (sign erased) or by lfcShrink(type='normal') (deprecated prior distorts the ranking). apeglm/ashr-shrunk results DROP the stat column - pull stat from the unshrunk results(dds) if ranking by it.

Run Preranked GSEA on GO

Goal: Find GO terms whose members shift coordinately up or down across the full ranking, with no significance cutoff.

Approach: Set the permutation seed for reproducibility, set eps=0 for exact tiny p-values, run gseGO, then map the leading-edge IDs back to symbols and read core_enrichment as the interpretable core.

set.seed(123)                              # fixes the multilevel Monte Carlo; any fixed seed
gse_go <- gseGO(geneList = gene_list, OrgDb = org.Hs.eg.db, keyType = 'ENTREZID',
                ont = 'BP', exponent = 1, minGSSize = 10, maxGSSize = 500,
                eps = 0, pvalueCutoff = 0.05, pAdjustMethod = 'BH',
                seed = TRUE, by = 'fgsea', verbose = FALSE)
gse_go <- setReadable(gse_go, OrgDb = org.Hs.eg.db, keyType = 'ENTREZID')

gseaResult columns: ID, Description, setSize, enrichmentScore (raw ES), NES, pvalue, p.adjust (BH), qvalue, rank (ES-peak position), leading_edge, core_enrichment (/-separated leading-edge IDs). Report p.adjust/qvalue, never raw pvalue and never an invented $FDR. NES sign = direction: positive = top of the ranking (up in the contrast), negative = bottom.

Run GSEA on KEGG, Reactome, or MSigDB

Goal: Apply the same preranked engine to a pathway database with the gene-ID type that database expects.

Approach: gseKEGG/gsePathway need ENTREZ-style IDs; choose the collection, keep the seed and eps=0, and note KEGG queries the live REST API (date-dependent) while Reactome and MSigDB are local.

set.seed(123)
gse_kegg <- gseKEGG(geneList = gene_list, organism = 'hsa', keyType = 'ncbi-geneid',   # Entrez-named vector; 'kegg' keyType is the prokaryote locus-tag path
                    minGSSize = 10, maxGSSize = 500, eps = 0,
                    pvalueCutoff = 0.05, seed = TRUE, verbose = FALSE)   # live KEGG API; pin the date

library(msigdbr)
h <- msigdbr(species = 'Homo sapiens', collection = 'H')    # 26.x: collection= (was category=); gs_collection (was gs_cat)
t2g <- h[, c('gs_name', 'ncbi_gene')]                       # 26.x Entrez column is ncbi_gene; older releases used entrez_gene
gse_h <- GSEA(geneList = gene_list, TERM2GENE = t2g, exponent = 1,
              minGSSize = 10, maxGSSize = 500, eps = 0,
              pvalueCutoff = 0.05, seed = TRUE, verbose = FALSE)

If the installed msigdbr still uses category=/entrez_gene, the old form works but warns - check ?msigdbr and names(h). ReactomePA's gsePathway(geneList, organism='human') reads the local reactome.db and also needs ENTREZ IDs.

Per-Sample Scores: ssGSEA and GSVA (Not a Contrast Test)

Goal: Convert an expression matrix into a gene-set-by-sample activity matrix to feed clustering, survival, or a classifier - there is no per-pathway p-value here.

Approach: GSVA >= 1.50 uses a PARAMETER-OBJECT API: build gsvaParam(...) or ssgseaParam(...) and pass it to gsva(). The old gsva(expr, gset.idx.list, method=) signature is defunct.

# GSVA >= 1.50 / Bioc 3.18 parameter-object API (older method= signature errors)
library(GSVA)
gsva_scores  <- gsva(gsvaParam(expr_matrix, gene_sets))      # unsupervised per-sample set scores
ssgsea_scores <- gsva(ssgseaParam(expr_matrix, gene_sets))   # ssGSEA via the same dispatch

Use GSEA (preranked or phenotype) for a CONTRAST and a pathway-level p-value; use ssGSEA/GSVA for a per-sample activity matrix for downstream modeling. GSVA is not installed in the reference environment - verify the installed signature with ?gsva before running.

Per-Method Failure Modes

Ranking by raw p-value

Trigger: gene_list <- -log10(de$pvalue) with no sign(). Mechanism: the magnitude is symmetric, so up- and down-regulated genes both land at the top. Symptom: NES signs are meaningless; "enriched" sets mix directions. Fix: rank by sign(log2FC) * -log10(pmax(p, 1e-300)), or use DESeq2 stat / limma t.

Preranked p-values treated as correlation-honest

Trigger: reporting FDR 0.001 from gseGO/fgsea on a co-regulated set. Mechanism: gene permutation assumes gene independence; correlated sets inflate the set-statistic variance, so p is too small. Symptom: "significant" pathways that are co-expression and do not replicate. Fix: state the permutation type; for type-I control with a design matrix use CAMERA (Wu & Smyth 2012).

Unsorted or duplicated geneList

Trigger: an un-sorted vector, or duplicate gene names after ID conversion. Mechanism: clusterProfiler assumes pre-sorting and uses names to map into sets; duplicates double-count a gene in the hit increments. Symptom: silently wrong ES, or an fgsea ties warning. Fix: sort(gl[!duplicated(names(gl))], decreasing=TRUE); prefer a continuous metric (Wald stat / moderated t rarely tie).

Bare log2FC ranking

Trigger: ranking by log2FoldChange from raw counts. Mechanism: a gene with 2 vs 8 counts shows a huge unstable LFC. Symptom: the leading edge is one or two low-count outliers, not a coordinated shift. Fix: rank by stat/t; if LFC is unavoidable use apeglm/ashr-shrunk LFC (never type='normal').

No set.seed

Trigger: running gseGO without fixing the seed. Mechanism: the multilevel Monte Carlo is stochastic. Symptom: p-values and the significant-set list drift across identical reruns. Fix: set.seed(123) AND seed=TRUE in the call.

Tiny leading edge believed

Trigger: trusting a high |NES| without inspecting core_enrichment. Mechanism: a 1-2 gene leading edge is outlier-driven, not a pathway shift; large sets reach high |NES| by chance. Symptom: an unreplicated headline pathway. Fix: FDR first, then leading-edge size/concentration, then NES for prioritization.

Stale nPerm / FDR<0.25 lore

Trigger: copying nPerm=10000 and "FDR < 0.25" from a Broad-desktop tutorial. Mechanism: nPerm was REMOVED at the fgsea/multilevel switch; clusterProfiler p.adjust is BH, not the Broad empirical-null FDR that 0.25 was calibrated for. Symptom: an argument error (nPerm) or a mis-transplanted threshold. Fix: drop nPerm, govern tiny-p with eps, treat p.adjust as BH and pick a defensible cutoff (often 0.05).

Stale GSVA call

Trigger: gsva(expr, gene_sets, method='ssgsea'). Mechanism: GSVA >= 1.50 dispatches on a parameter object's class. Symptom: the old signature errors. Fix: gsva(gsvaParam(expr, gene_sets)) / ssgseaParam(...).

Quantitative Thresholds

| Threshold | Source | Rationale | |-----------|--------|-----------| | exponent = 1 | Subramanian 2005 PNAS 102:15545 | weights each hit by |statistic|; p=0 is the unweighted KS that flags middle-clustered sets with no strong genes | | minGSSize = 10 | clusterProfiler default | drops tiny sets that overfit on one outlier | | maxGSSize = 500 | clusterProfiler default | drops overly broad sets that always 'enrich' | | eps = 0 (default 1e-10) | clusterProfiler / fgsea | replaces nPerm: eps=0 resolves exact tiny p-values; 1e-10 is the default floor | | pAdjustMethod = 'BH' | clusterProfiler default | Benjamini-Hochberg FDR; NOT the Broad empirical-null FDR, so do not reflex to 0.25 | | pvalueCutoff = 0.05 | clusterProfiler default | filters on p.adjust by default; defensible BH cutoff | | ~>=7 samples/group | Broad GSEA docs | minimum for a non-degenerate phenotype-permutation null | | set.seed(123) | reproducibility | any fixed seed; the point is to fix the multilevel Monte Carlo |

Common Errors

| Error / symptom | Cause | Solution | |-----------------|-------|----------| | Error about names / wrong ES | geneList not named or not sorted decreasing | sort(setNames(v, ids), decreasing=TRUE) | | --> No gene can be mapped | wrong keyType/OrgDb, or non-ENTREZ IDs | bitr to the expected ID type first | | gseKEGG returns 0 terms | ENSEMBL/SYMBOL passed, wrong organism code, or KEGG API down | convert to kegg-id/ENTREZ; check the organism code; retry (live API) | | Different results each run | no set.seed, or live KEGG DB changed | fix the seed; pin the KEGG run date | | nPerm argument error | copied from a pre-4.0 tutorial | remove nPerm; use eps | | GSVA method= error | pre-1.50 signature | gsva(gsvaParam(expr, sets)) | | core_enrichment is NA / all-ID | setReadable not applied | setReadable(gse, OrgDb, keyType='ENTREZID') |

References

• Mootha VK, Lindgren CM, Eriksson KF, et al. 2003. PGC-1alpha-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat Genet 34:267-273.
• Subramanian A, Tamayo P, Mootha VK, et al. 2005. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. PNAS 102:15545-15550.
• Korotkevich G, Sukhov V, Budin N, et al. Fast gene set enrichment analysis. bioRxiv 060012 (preprint). DOI 10.1101/060012.
• Wu D, Smyth GK. 2012. Camera: a competitive gene set test accounting for inter-gene correlation. Nucleic Acids Res 40:e133.
• Wu D, Lim E, Vaillant F, et al. 2010. ROAST: rotation gene set tests for complex microarray experiments. Bioinformatics 26:2176-2182.
• Barbie DA, Tamayo P, Boehm JS, et al. 2009. Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature 462:108-112.
• Hanzelmann S, Castelo R, Guinney J. 2013. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 14:7.
• Reimand J, Isserlin R, Voisin V, et al. 2019. Pathway enrichment analysis and visualization of omics data using g:Profiler, GSEA, Cytoscape and EnrichmentMap. Nat Protoc 14:482-517.
• Wijesooriya K, Jadaan SA, Perera KL, et al. 2022. Urgent need for consistent standards in functional enrichment analysis. PLoS Comput Biol 18:e1009935.

Related Skills

• go-enrichment - Gene-list ORA alternative when no ranking exists
• kegg-pathways - KEGG pathway and module enrichment and ID conventions
• reactome-pathways - Reactome curated-pathway enrichment (local DB)
• wikipathways - WikiPathways community-pathway enrichment
• enrichment-visualization - gseaplot2, ridgeplot, and dotplot of GSEA results
• differential-expression/de-results - Source of the ranking statistic and the padj column conventions
• workflows/expression-to-pathways - End-to-end DE-to-enrichment pipeline