🧪 Metabolic Pathway Enrichment Analysis

Untargeted Metabolomics — Capturing the Biofluid Chemical Landscape

Pathway enrichment analysis begins not with a hypothesis but with an unbiased chemical survey. Untargeted LC-MS/MS metabolomics profiles thousands of small molecules simultaneously — amino acids, lipids, organic acids, nucleotides, xenobiotics — without preselecting which ones matter. The quality of every downstream enrichment result depends entirely on how carefully this acquisition step is controlled, since batch effects and ion suppression introduced here propagate through statistics, mapping, and biological interpretation alike.

• 220: Cohort size (120 cases / 100 controls)
• QTOF: Instrument (Agilent 6546, orthogonal accel.)
• HILIC + C18: Chromatography (covers polar + lipophilic space)
• <3 ppm: Mass accuracy (enables putative ID by m/z)

Sample handling and extraction protocol

Plasma metabolomics is exquisitely sensitive to pre-analytical variation, so protocol standardization is the first line of defense against confounding:

• Collection: blood drawn into EDTA tubes, centrifuged at 1,300×g for 10 min at 4°C within 30 minutes of draw; plasma aliquoted and snap-frozen at −80°C • Extraction: cold methanol (−20°C) added at 1:4 plasma:solvent ratio to simultaneously precipitate proteins and quench enzymatic activity; vortexed, incubated 10 min at −20°C, centrifuged 15,000×g • QC pooling: a pooled QC sample (equal aliquot from every study sample) is injected every 8–10 runs to monitor instrument drift and enable signal-drift correction • Randomization: case/control samples are injected in randomized, blocked order to prevent batch effects from being confounded with disease status

Extraction efficiency and matrix effects are evaluated with spiked internal standards (e.g., deuterated amino acids, ¹³C-labeled fatty acids) added before extraction — recovery outside 70–120% for any standard flags a batch for re-extraction.

Chromatographic separation and dual-polarity mass spectrometry

No single chromatographic mode captures the full metabolome, so paired columns are run to maximize coverage:

• HILIC (hydrophilic interaction, e.g., Waters BEH Amide): retains and separates highly polar species — amino acids, sugars, nucleotides, organic acids — that elute unretained on reversed-phase columns • C18 reversed-phase: separates lipophilic species — phospholipids, sphingolipids, acylcarnitines, bile acids, steroids • Gradient: 18-minute run per column, 0.4 mL/min, mobile phases of water/acetonitrile with 0.1% formic acid (positive mode) or ammonium acetate (negative mode) • Dual ESI polarity: positive mode ionizes basic/neutral species efficiently (amines, many lipids); negative mode favors acidic species (fatty acids, nucleotides, most sugars) — running both roughly doubles metabolome coverage

Each sample therefore generates four raw data files (HILIC+/HILIC−/C18+/C18−), each containing on the order of 2,000 chromatographic features, for a combined ~8,400 raw features per subject before any deduplication or annotation.

Peak Picking, Alignment, and Differential Abundance Testing

Raw chromatographic data is dominated by noise: isotopes, in-source fragments, adducts, and instrument artifacts vastly outnumber true unique metabolites. Computational processing collapses this raw feature table into a defensible list of annotated, quantified molecules, and rigorous statistical testing then narrows that list to the handful of species that actually differ between biological groups — the input list that pathway enrichment will operate on.

• centWave: Peak-picking algorithm (XCMS3, wavelet-based)
• obiwarp: Alignment method (nonlinear RT correction)
• Mann–Whitney U: Test statistic (non-parametric, per-metabolite)
• Benjamini–Hochberg: Multiple testing (controls FDR across 1,240 tests)

From raw signal to annotated metabolite table

Untargeted LC-MS data processing is a multi-stage computational pipeline, typically executed in XCMS3, MZmine 3, or MS-DIAL:

1. Peak picking (centWave): detects chromatographic peaks in the m/z–retention time plane using continuous wavelet transforms tuned to expected peak width (typically 5–20 seconds) and mass accuracy tolerance (5–15 ppm) 2. Retention time alignment (obiwarp): corrects run-to-run RT drift using a nonlinear warping function anchored to the pooled QC injections, ensuring the same metabolite is matched across all 220 samples 3. Grouping/correspondence: features with matching aligned RT and m/z across samples are grouped into a single "feature" row in the final data matrix 4. CAMERA/annotation: groups isotope clusters (¹³C satellites), common adducts ([M+H]⁺, [M+Na]⁺, [M−H]⁻), and in-source fragments belonging to the same parent compound, collapsing redundant signals 5. Putative identification: surviving features are matched against spectral libraries (mzCloud, METLIN, HMDB) by accurate mass (±5 ppm) and, where MS/MS was collected, fragmentation spectral match score

After this pipeline, ~8,400 raw features per polarity/column combination collapse to 1,240 confidently annotated, non-redundant metabolites suitable for statistical testing.

Statistical filtering for differential abundance

Each of the 1,240 metabolites is tested independently for a case-vs-control abundance difference:

• Normalization: total ion current (TIC) normalization or probabilistic quotient normalization (PQN) corrects for dilution differences between plasma samples; log2 transformation stabilizes variance • Univariate test: Mann–Whitney U (non-parametric, robust to the non-normal distributions typical of MS intensity data) compares each metabolite between groups; parametric t-tests are used when normality holds • Effect size: log2 fold-change (log2FC) between group medians quantifies the magnitude, independent of statistical significance • Multiple testing correction: with 1,240 simultaneous tests, uncorrected p<0.05 would yield ~62 false positives by chance alone — Benjamini–Hochberg FDR correction controls the expected proportion of false discoveries among all "significant" calls • Final hit list: metabolites passing FDR<0.05 AND |log2FC|>0.5 (roughly 1.4-fold change) are retained — in this cohort, 187 metabolites meet both criteria

This 187-metabolite hit list — not the full 1,240-metabolite background — is the direct input to pathway mapping and enrichment testing in the next stages.

Projecting Hit Metabolites onto KEGG and Reactome Pathway Maps

A list of 187 altered metabolite names or masses carries little biological meaning on its own. Pathway mapping translates each hit into a formal compound identifier and locates it within curated metabolic network databases — KEGG and Reactome — that encode which enzymes, reactions, and pathways each molecule participates in. This step converts an unstructured hit list into a structured biological network that statistics can then be run over.

• ~350: KEGG human pathways (metabolic + signaling total)
• 84: Pathways queried here (metabolic pathway subset)
• 3: Compound ID systems (KEGG C-number, HMDB, ChEBI)
• 74%: Mapping success rate (139 of 187 hits placed)

Compound identifier cross-referencing

The central obstacle in pathway mapping is that mass spectrometry identifies molecules by mass and fragmentation pattern, while pathway databases index molecules by curated compound identifiers — these vocabularies must be reconciled:

• KEGG Compound (C-numbers, e.g., C00025 = L-glutamate): the reference vocabulary for KEGG pathway maps; each metabolite hit is matched by exact name, InChIKey, or monoisotopic mass • HMDB (Human Metabolome Database) accessions: provide cross-links to KEGG, ChEBI, and PubChem, and are the most complete resource for human-specific metabolites including many lipids absent from KEGG • ChEBI (Chemical Entities of Biological Interest): used primarily for Reactome pathway cross-referencing, with a more granular ontology of chemical classes • Tools: MetaboAnalyst 6.0 compound name-matching utility and the KEGG API automate this cross-referencing at scale, but ambiguous or trivial metabolite names (e.g., isomeric sugars, unresolved lipid species written only as sum-composition like "PC 34:2") often require manual curation

Of 187 hits, 139 (74%) resolve unambiguously to a mapped compound ID; the remaining 48 are either genuinely absent from curated pathway databases (many complex lipids and xenobiotic metabolites) or are annotated only putatively and excluded from pathway-level testing to avoid introducing false mappings.

Pathway database architecture — KEGG vs. Reactome

KEGG and Reactome encode metabolic knowledge differently, and combining both improves coverage and topological resolution:

• KEGG pathway maps: hand-curated reaction diagrams (e.g., hsa00250 Alanine/aspartate/glutamate metabolism) representing metabolites as nodes and enzymatic reactions as edges; pathways are relatively coarse-grained (median ~40 compounds per pathway) and organized by classical biochemical pathway names • Reactome: a more granular, reaction-centric database representing individual enzymatic steps, complexes, and regulatory events; a single KEGG pathway may correspond to several more specific Reactome sub-pathways, giving finer resolution for topology-based scoring • Pathway size distribution matters statistically: very large pathways (e.g., "metabolic pathways," >1,000 compounds) are typically excluded or down-weighted, since almost any hit list will show nominal overlap with them by chance; very small pathways (<3 compounds) are also excluded because enrichment statistics are unstable with so few members • In this analysis, the queried set of 84 pathways was restricted to those with 3–200 member compounds annotated in the human metabolome, and at least one member among the 139 mapped hits

Over-Representation Analysis and Topology-Weighted Scoring

With hit metabolites mapped onto pathway networks, the core statistical question becomes: are more hits falling into a given pathway than expected by chance, and do those hits cluster at biologically influential positions within the pathway's reaction network? Two complementary strategies answer this — classical over-representation analysis (ORA) and topology-aware algorithms like mummichog that weight hits by their network centrality.

• Hypergeometric: Primary test (ORA, one-sided)
• Mummichog v2: Topology algorithm (Li et al., PLOS Comp Bio 2013)
• 1,000 permutations: Empirical p-value (randomized hit sets)
• 9 / 84: Significant pathways (FDR<0.01)

Hypergeometric over-representation testing

ORA asks a simple combinatorial question for each pathway: given a background of N total mapped compounds, K of which belong to this pathway, what is the probability of drawing k or more pathway members by chance among the n hits actually observed?

P(X≥k) = Σ [C(K,i)·C(N−K,n−i)] / C(N,n), summed for i=k…min(n,K)

Worked example — arginine and proline metabolism (K=32 compounds in a background of N=1,240 mapped metabolites): among the n=139 mapped hits, k=11 belong to this pathway. The hypergeometric test gives p=3.1×10⁻⁵ — far more overlap than the ~3.6 hits expected by chance alone (139×32/1240).

ORA limitations that motivate topology-aware alternatives: • Treats every pathway member as equally important, ignoring network position — a hit on a rate-limiting bottleneck enzyme carries more biological weight than one on a peripheral branch • Requires a hard hit/non-hit cutoff (the FDR<0.05, |log2FC|>0.5 threshold from Stage 2), discarding continuous evidence from metabolites just below threshold • Pathway overlap (shared compounds between related pathways, e.g., multiple amino acid pathways sharing glutamate) inflates redundant "hits," which is why FDR correction is applied across all 84 tested pathways, not per-test

Mummichog — topology and permutation-based significance

Mummichog (Li et al., PLOS Computational Biology 2013) was designed specifically for untargeted metabolomics, where compound identification is uncertain, by testing pathway enrichment using putative annotations directly from m/z values combined with network topology:

• Local network activity score: for each pathway, a composite score combines (a) the fraction of pathway members among the hit list and (b) the betweenness centrality of hit-metabolite nodes within the full metabolic reaction network — hits at network chokepoints score higher than hits at dead-end branches • Permutation testing: because centrality distributions are not analytically tractable, empirical null distributions are built by repeatedly (1,000×) drawing random hit sets of the same size from the full annotated metabolome and recomputing the activity score — the fraction of permuted scores exceeding the observed score gives the empirical p-value • Pathway impact score: for the top pathway in this analysis (arginine biosynthesis), the topology-weighted impact score reaches 0.42 on a 0–1 scale, reflecting that hit metabolites sit on several branch points controlling flux through the pathway, not merely at peripheral positions • Combined with ORA and corrected across 84 tested pathways by Benjamini–Hochberg FDR, 9 pathways remain significant at FDR<0.01 — led by arginine/proline metabolism, glycine/serine/threonine metabolism, and the tryptophan-kynurenine pathway

From Enriched Pathways to a Validated Clinical Biomarker Panel

Statistical enrichment identifies which pathways are perturbed, but converting that into clinical utility requires a further translational step: nominating a compact biomarker panel from the enriched pathways, and confirming its discriminative performance with targeted, quantitative assays in an independent patient cohort. This final stage is where computational pathway biology meets diagnostic validation.

• 5: Nominated biomarkers (from top-3 enriched pathways)
• n=80: Validation cohort (independent of discovery set)
• MRM-MS: Targeted assay (triple-quadrupole, absolute quant.)
• 0.89: Panel AUC (95% CI 0.83–0.94)

From pathway enrichment to biomarker panel selection

The 9 significantly enriched pathways collectively implicate dozens of metabolites, far too many for a practical clinical assay — panel selection narrows this to a small, mechanistically coherent, and analytically robust set:

• Mechanistic coherence: metabolites are prioritized if they sit at pathway bottlenecks flagged by the mummichog impact score, since these are more likely to reflect causal perturbation of enzyme activity rather than passive downstream drift • Analytical robustness: candidates must have a commercially available stable-isotope-labeled internal standard for absolute quantification, and a validated MRM transition free of isobaric interference • Effect size and reproducibility: candidates are ranked by consistency of log2FC and significance across the full discovery cohort, not just the pooled test • Final panel: kynurenine (tryptophan-kynurenine pathway, elevated), propionylcarnitine/C3-acylcarnitine (branched-chain amino acid catabolism, elevated), lysophosphatidylcholine LPC 18:2 (glycerophospholipid metabolism, depleted), glutamate (arginine/proline metabolism, elevated), and 2-hydroxybutyrate (glutathione/cysteine metabolism, elevated) — five metabolites spanning three of the top enriched pathways

Independent cohort validation with targeted mass spectrometry

Untargeted discovery data is semi-quantitative and prone to batch effects, so clinical translation requires re-measuring the nominated panel with a targeted, absolutely quantitative method in new patients:

• Assay: multiple reaction monitoring (MRM) on a triple-quadrupole LC-MS/MS instrument, with two transitions per compound (one quantifier, one qualifier) and a stable-isotope-labeled internal standard for each analyte • Cohort: 80 new subjects (independent of the 220-subject discovery cohort), recruited from a separate clinical site to test generalizability and reduce batch/site confounding • Model: logistic regression combining the five standardized metabolite concentrations, with performance assessed by 10-fold cross-validation and reported via ROC analysis • Result: the 5-metabolite panel achieves AUC 0.89 (95% CI 0.83–0.94) for case/control classification, outperforming any single metabolite (best individual AUC: kynurenine at 0.76) and approaching the performance of the full 187-metabolite discovery signature (AUC 0.91) with a fraction of the analytical cost

This pattern — pathway enrichment nominating a small mechanistically grounded panel that then matches most of the discriminative power of the full untargeted signature — is the typical and desired outcome of a well-executed enrichment workflow.

The five-metabolite panel derived from pathway enrichment (AUC 0.89) recovers 98% of the discriminative power of the full 187-metabolite untargeted signature (AUC 0.91) at roughly 1/40th the analytical cost per sample. This demonstrates the core value proposition of pathway enrichment analysis: it is not merely a bioinformatic annotation exercise, but a data-driven method for compressing thousands of untargeted features into a small, mechanistically interpretable, clinically deployable panel.