Bio Pathway Reactome

Version Compatibility

Reference examples tested with: ReactomePA 1.54+, reactome.db 1.95+, clusterProfiler 4.18+.

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.

reactome.db is a LOCAL Bioconductor annotation package pinned to the Bioconductor release, so enrichPathway/gsePathway are reproducible offline given the package version - unlike KEGG and WikiPathways, which query a live database. The reactome.org web AnalysisService tracks the current quarterly Reactome release and can therefore disagree with a local ReactomePA run (see Failure Modes).

Reactome Pathway Enrichment

"Which curated Reactome pathways does my gene list over-represent?" -> Test each Reactome pathway for over-representation against a measured background, then deduplicate the hierarchy - because a Reactome result is one signal projected onto a tree of nested reactions, not a list of independent findings.

• R: enrichPathway(gene_entrez, organism='human', universe=measured_entrez, readable=TRUE)

Scope: ORA (enrichPathway) and GSEA (gsePathway) over Reactome reaction-rolled-to-pathway gene sets, the ENTREZ-only constraint, hierarchy deduplication, the human-curated-only species caveat, the local viewPathway reaction-network plot, and the ReactomeGSA comparative pointer. The hypergeometric test and background-universe theory -> go-enrichment. The GSEA running-sum engine and ranking metric -> gsea. The DE list / ranking statistic -> differential-expression/de-results. Dotplot/emapplot/cnetplot/gseaplot2 -> enrichment-visualization.

The Single Most Important Modern Insight -- Reactome's Atomic Unit Is the Reaction and Pathways Are Nested Containers, So a Result Is One Signal Projected Onto a Tree, Not a List

Reactome is not a collection of pathway maps like KEGG. Its atomic unit is the ReactionlikeEvent - a single typed molecular transformation (binding, catalysis, transport, modification) with a PubMed citation - and pathways are containers that group reactions and sub-pathways into a deep event hierarchy (TopLevelPathway -> Pathway -> sub-Pathway -> Reaction). Two consequences define every decision in this skill:

1. Granularity is the reason to choose Reactome AND the multiple-testing tax. A Reactome hit is finer than a KEGG map - a specific, peer-reviewed, literature-grounded reaction - which is why it is worth using. But finer means MORE gene sets (many tiny leaf pathways plus a few huge top-level ones), so the multiple-testing burden is heavier and minGSSize/maxGSSize matter more than for GO.
2. Nesting double-counts the signal. A gene annotated to one leaf reaction is a member of that pathway AND every ancestor, so a parent and child enrich on the SAME genes. A live cell-cycle gene list returns "Cell Cycle Checkpoints" (parent), "G2/M Checkpoints", and "G1/S Transition" (children) stacked at the top - one signal, three "independent" small p-values. A Reactome table is therefore read by reasoning about the hierarchy (report the deepest significant node, ancestors as context), not by sorting on p.adjust.

Second load-bearing fact: only human is curated; every non-human pathway is orthology-inferred from the human reactions, not independently curated. Mouse projection is ~81% complete, exotic species far less, so a "Reactome mouse pathway" is a hypothesis from orthology that inherits human-curation gaps and misses mouse-specific biology. The honest output is "the deepest curated human pathway nodes my list over-represents, deduplicated against their ancestors, against a background of genes I actually measured."

Tool Taxonomy

Decision Tree by Scenario

ReactomePA's organism accepts exactly seven values: human, rat, mouse, celegans, yeast, zebrafish, fly. This is a reactome.db mapping ceiling, NOT a Reactome ceiling - the database projects to ~14-20 species and the web AnalysisService covers them; do not conflate the two.

Over-Representation Analysis (enrichPathway)

Goal: Find Reactome pathways over-represented in a significant-gene list, against the genes actually measured.

Approach: Convert the gene list to ENTREZ (mandatory - no keyType argument), pass the measured background as universe, run enrichPathway, then read the hierarchy rather than the raw row order.

library(ReactomePA)
library(org.Hs.eg.db)
library(clusterProfiler)   # bitr

sig_entrez <- bitr(sig_symbols, fromType='SYMBOL', toType='ENTREZID', OrgDb=org.Hs.eg.db)$ENTREZID
universe   <- bitr(all_tested_symbols, fromType='SYMBOL', toType='ENTREZID', OrgDb=org.Hs.eg.db)$ENTREZID

ora <- enrichPathway(gene=sig_entrez, organism='human', universe=universe,
                     pvalueCutoff=0.05, qvalueCutoff=0.2,
                     minGSSize=10, maxGSSize=500, readable=TRUE)
# enrichResult columns: ID Description GeneRatio BgRatio RichFactor FoldEnrichment zScore
#                       pvalue p.adjust qvalue geneID Count
# FoldEnrichment and RichFactor are columns - read effect size there, do not compute GeneRatio/BgRatio by hand.

Without universe, the background is all ~11,200 Reactome-annotated ENTREZ genes (the live BgRatio denominator), not the genes measured, which over-states significance. readable=TRUE maps the geneID column back to symbols. Reading the result means deduplicating the hierarchy: identify the deepest significant pathway for each signal and note its ancestors as context, not as separate hits. ReactomePA has no simplify() equivalent (unlike GO's DAG), so this is a manual judgment call; the visual collapse (treeplot, emapplot) is owned by enrichment-visualization.

GSEA (gsePathway)

Goal: Find Reactome pathways whose genes shift coordinately across the full ranking, without a significance cutoff.

Approach: Build a named numeric vector sorted decreasing by the ranking statistic with ENTREZ names, set a seed for permutation reproducibility, then run gsePathway and read the leading edge.

gene_list <- de$stat                       # any per-gene statistic: t-stat, signed -log10 p, shrunken log2FC
names(gene_list) <- de$entrez              # names MUST be ENTREZ
gene_list <- sort(gene_list, decreasing=TRUE)

set.seed(123)                              # gsePathway permutes; fix the seed so p-values reproduce
gse <- gsePathway(geneList=gene_list, organism='human',
                  pvalueCutoff=0.05, pAdjustMethod='BH', verbose=FALSE)
# gseaResult columns: ID Description setSize enrichmentScore NES pvalue p.adjust qvalue rank leading_edge core_enrichment

GSEA uses the whole ranking, so the universe/background pitfall of ORA does not apply - but the hierarchy double-counting STILL does: a parent and child both score on the same leading-edge genes. The ranking metric IS the experiment (the same genes ranked differently give different leading edges); the metric choice is owned by gsea.

viewPathway - the Reactome reaction-network plot

Goal: Draw the reactions and physical entities of ONE enriched pathway, optionally colored by fold change.

Approach: Pass the pathway NAME (the Description, NOT the R-HSA id) to viewPathway; it renders a LOCAL ggraph reaction network in the R graphics device. To open the actual web diagram, build the PathwayBrowser URL from the R-HSA id and browseURL it.

top_name <- ora@result$Description[1]                       # the NAME, not $ID
viewPathway(top_name, organism='human', readable=TRUE, foldChange=gene_list)

# the interactive web diagram needs the R-HSA id, NOT viewPathway:
browseURL(paste0('https://reactome.org/PathwayBrowser/#/', ora@result$ID[1]))

viewPathway's keyType controls the ID type of the foldChange names only (so a SYMBOL-named fold-change vector works with keyType='SYMBOL'); the ENTREZ-only constraint of enrichPathway does not extend to it. Route generic dotplot/emapplot/cnetplot/gseaplot2 to enrichment-visualization; viewPathway is owned here because it is Reactome-data-structure-specific.

ReactomeGSA - comparative and multi-omics

Goal: Compare pathway activity BETWEEN conditions, across omics layers, or across single-cell clusters - which ReactomePA's single-list model cannot do.

Approach: ReactomeGSA is a separate Bioconductor client to Reactome's hosted AnalysisService; build a request, add datasets, and send it to the server (network required; results track the server's release, not a local reactome.db).

library(ReactomeGSA)
req <- ReactomeAnalysisRequest(method='Camera')                  # or 'ssGSEA', 'PADOG'
req <- add_dataset(req, expression_values=expr_matrix, name='RNAseq',
                   type='rnaseq_counts', comparison_factor='condition',
                   comparison_group_1='A', comparison_group_2='B', sample_data=meta)
res <- perform_reactome_analysis(req)                            # sends to the server
pw  <- pathways(res)                                             # combined pathway table
sc  <- analyse_sc_clusters(seurat_obj, use_interactors=FALSE)    # per-cluster ssGSEA

Use ReactomePA for "is this one list over-represented / coordinately changed"; use ReactomeGSA for "which pathways DIFFER between conditions / omics / clusters".