Introduction

Advances in metabolomics allow hundreds of metabolites to be measured in a single experiment, but linking those features to biology still depends on pathway enrichment, set enrichment, and network analysis. Existing tools often chain multiple steps, rely on web services, or are hard to install reproducibly. enrichmet is an R/Bioconductor package for local, streamlined pathway analysis of metabolomics data. A single call to enrichmet() can run:

over-representation analysis (Fisher’s exact test);
metabolite set enrichment analysis (MetSEA) via fgsea;
pathway–metabolite network topology (centrality via igraph);
pathway–metabolite and Reactome-based interaction visualisations.

The workflow is modular: internal steps are available as separate functions with examples, so users can run the full pipeline or customise individual stages. Enrichment runs locally (no required web UI for the statistical workflow), which suits high-throughput and iterative analyses.

Pathway and compound annotations are obtained at runtime from public APIs (KEGG REST, Reactome mappings, optional LION/BioPortal) and cached with BiocFileCache. The package does not redistribute third-party database dumps (KEGG, Reactome, or LION). Equivalent Pathway / Metabolites tables may be supplied by the user at runtime when preferred.

enrichmet() returns result tables (pathway enrichment, MetSEA, centrality) and, optionally, plots including pathway enrichment, impact, MetSEA, relative betweenness (RBC), pathway–metabolite network, heatmap, membership, and Reactome interaction network. If differential analysis is run with run_de(), a volcano plot can also be included.

Installation

if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}
BiocManager::install("enrichmet")

library(enrichmet)

Input requirements

enrichmet() needs a query set of metabolites and, for MetSEA or a measured background, summary statistics.

Query metabolites

Provide one of:

inputMetabolites — character vector of KEGG IDs and/or names (see Section 3.3); or
da_results — output of run_de(), from which the query set is derived.

`input_type`	Behaviour
`"kegg"` (default)	Entries treated as KEGG IDs (`C#####`)
`"name"`	Entries mapped to KEGG IDs via `kegg_lookup`
`"mixed"`	Valid KEGG IDs kept; remaining entries mapped as names

For "name" or "mixed", supply kegg_lookup (e.g. from fetch_kegg_compound_lookup()).

Summary statistics (optional but required for MetSEA)

A minimum run needs only the query set above. For MetSEA (and for building a measured-metabolite background), also provide a table in the format of example_data (Section 3.3) with:

met_id — feature identifier (may embed KEGG IDs);
pval or padj — statistical significance;
log2fc — log2 fold change.

MetSEA ranks metabolites by effect size and statistical evidence, so enrichment uses the ordered list rather than a binary significant / non-significant split.

Pathway–metabolite map

enrichmet requires a data frame with columns Pathway and Metabolites (comma-separated compound identifiers). This table is the annotation background: it defines which metabolites belong to which pathways and serves as the reference universe for Fisher’s exact test, so pathway over-representation can be assessed against known pathway membership.

For KEGG, mappings are not shipped with the package. They are retrieved at runtime from the KEGG REST API (https://rest.kegg.jp/) when the user calls fetch_kegg_pathway_metabolites(). Responses are cached locally with BiocFileCache so repeated calls in the same environment do not re-download unnecessarily.

The organism argument is a KEGG organism code. For example, human (hsa):

# Fetches human pathway list and compound–pathway links via the KEGG REST API
PathwayVsMetabolites <- fetch_kegg_pathway_metabolites(organism = "hsa")

Case study

This case study uses human KEGG pathway maps and compound names retrieved via the KEGG REST API, example summary statistics for demonstration, and Reactome reaction mappings built at runtime.

inputMetabolites <- c(
    "C00209", "C00249", "C06424", "C01727", "C04025", "C00366", "C00385",
    "C08261", "C00079", "C00506", "C00093", "C02979", "C00106", "C01384",
    "C00219", "C00242", "C19806", "C00105", "C00606", "C00144", "C01216",
    "C06429", "C00475", "C00386", "C01046", "C02180", "C00082", "C00299",
    "C00051", "C16513", "C00387", "C00022", "C00026", "C01062", "C00418",
    "C00295", "C01762", "C00525", "C00074", "C00491", "C02067", "C00015",
    "C14829", "C00519", "C06428", "C05472", "C00328", "C16358", "C16357",
    "C16353", "C02632", "C00246", "C05842", "C00785", "C03739", "C00212",
    "C01181", "C00140", "C05021", "C02862", "C02721", "C01104", "C10172",
    "C00380", "C01262", "C05122", "C03736", "C00599", "C00446", "C00275",
    "C00103", "C00881", "C15587", "C00361", "C00429", "C01042", "C11439",
    "C05637", "C00127", "C00035", "C01909", "C00836", "C00847", "C04294",
    "C00043", "C01481", "C00301", "C05635", "C02140", "C05488", "C00570",
    "C01551", "C00092", "C03139"
)

PathwayVsMetabolites <- fetch_kegg_pathway_metabolites(organism = "hsa")
#> Filter applied: excluded 280 pathways. Remaining: 92
kegg_lookup <- fetch_kegg_compound_lookup()

example_path <- system.file(
    "extdata", "summary_stat.csv",
    package = "enrichmet"
)

if (example_path == "") {
    stop("Example file 'summary_stat.csv' not found in inst/extdata/")
}

example_data <- read.csv(example_path, stringsAsFactors = FALSE)

kegg_ids <- extract_kegg_ids(PathwayVsMetabolites)
#> Pathway universe contains 2960 unique KEGG compound IDs
reactome_df <- fetch_kegg_reactome(
    kegg_ids = kegg_ids,
    species = "Homo sapiens"
)

results <- enrichmet(
    inputMetabolites = inputMetabolites,
    PathwayVsMetabolites = PathwayVsMetabolites,
    example_data = example_data,
    kegg_lookup = kegg_lookup,
    reactome_df = reactome_df,
    top_n = 20,
    p_value_cutoff = 0.5,
    network_top_n = 10,
    heatmap_top_n = 20,
    membership_top_n = 20,
    min_pathway_occurrence = 2,
    min_metabolite_occurrence = 1,
    input_type = "kegg"
)
#> Using 94 metabolites from inputMetabolites character vector
#> Processing complex KEGG IDs (splitting by |)...
#> Split 94 complex IDs into 94 unique KEGG IDs
#> Running pathway enrichment analysis...
#> Enrichment analysis completed: 60 pathways tested, 60 pathways passed filtering (p <= 1.000000)
#> Running GSEA analysis...
#> Extracting KEGG IDs from 'met_id' column...
#> Successfully extracted KEGG IDs for 398 metabolites
#> Sample extracted KEGG IDs: C00245, C01879, C00025, C00209, C00158, C00049
#> Prepared 243 metabolites with KEGG IDs for GSEA
#> Created rankings for 243 KEGG metabolites using 'signed_pval' method
#> Ranking range: -10.812 to 8.832
#> Preparing GMT data from 84 pathways
#> GMT preparation summary:
#>   Total pathways processed: 84
#>   Pathways kept (size >= 5): 40
#>   Pathways removed (size < 5): 44
#>   Pathway size distribution:
#>     Min: 5
#>     Median: 7
#>     Max: 21
#>     Mean: 8.9
#> Testing 40 pathways with GSEA
#> GSEA completed: 40 pathways tested, 1 significant at FDR < 0.05
#> Top pathway: 'Aminoacyl-tRNA biosynthesis' (NES = -1.89, FDR = 0.048)
#> Applied KEGG pathway name mapping using 'kegg_id' and 'name' columns
#> Running centrality analysis...
#> Generating metabolite-pathway network visualization...
#> Using top 10 metabolites by centrality
#> Generating enrichment heatmap...
#> Selected top 20 metabolites by enrichment significance for heatmap
#> Generating pathway membership plot...
#> === MEMBERSHIP PLOT DEBUG ===
#> Input metabolites: 94
#> Top_n parameter: 20
#> Initial matched metabolites: 72
#> Filtering to top 20 metabolites by centrality...
#> Total metabolites with centrality: 2961
#> Top metabolites selected: 20
#> Sample top metabolites: C00022, C00026, C00103, C00082, C00249, C00301
#> After centrality filtering - unique metabolites: 20
#> Matrix before occurrence filtering: 20 x 51
#> Pathway counts range: 1 to 4
#> Metabolite counts range: 2 to 24
#> Pathways to keep: 30/51
#> Metabolites to keep: 20/20
#> Final matrix dimensions: 20 x 30
#> === END DEBUG ===
#> Generating Reactome interaction network...
#> Extracting KEGG IDs from input metabolites...
#> Successfully extracted 94 KEGG IDs from 94 input metabolites
#> Sample extracted KEGG IDs: C00209, C00249, C06424, C01727, C04025, C00366
#> Found 25 metabolites with 292 associated Reactome reactions
#> Building metabolite co-occurrence edges from shared reactions...
#> Found 31 valid metabolite pairs sharing >= 1 reaction(s)
#> Applied KEGG pathway name mapping
#> Display names sample: Pyruvate, 2-Oxoglutarate, Uracil, Adenosine, Guanine, Butanoic acid
#> Creating graph with 15 vertices
#> Vertex attributes: name, display_name
#> Graph vertex attributes: name, display_name
#> Sample vertex display_names: Pyruvate, 2-Oxoglutarate, Uracil, Adenosine, Guanine, Butanoic acid
#> Using layout: gem (spacing score: 425.94)

Background-corrected enrichment

Fisher’s exact test can be restricted to metabolites that were measured (testable) in the experiment, reducing bias from database compounds that were never observed.

measured_background <- extract_measured_kegg_ids(example_data)

results_bg <- enrichmet(
    inputMetabolites = inputMetabolites,
    PathwayVsMetabolites = PathwayVsMetabolites,
    example_data = example_data,
    kegg_lookup = kegg_lookup,
    reactome_df = reactome_df,
    top_n = 20,
    p_value_cutoff = 1,
    network_top_n = 10,
    heatmap_top_n = 20,
    membership_top_n = 20,
    min_pathway_occurrence = 2,
    min_metabolite_occurrence = 1,
    backgroundMetabolites = measured_background,
    input_type = "kegg"
)
#> Using 94 metabolites from inputMetabolites character vector
#> Processing complex KEGG IDs (splitting by |)...
#> Split 94 complex IDs into 94 unique KEGG IDs
#> Running pathway enrichment analysis (background-corrected)...
#> Background correction applied: universe restricted from 2960 (all pathway-annotated metabolites) to 214 (measured metabolites also present in the pathway database)
#> Enrichment analysis completed: 60 pathways tested, 60 pathways passed filtering (p <= 1.000000)
#> Running GSEA analysis...
#> Extracting KEGG IDs from 'met_id' column...
#> Successfully extracted KEGG IDs for 398 metabolites
#> Sample extracted KEGG IDs: C00245, C01879, C00025, C00209, C00158, C00049
#> Prepared 243 metabolites with KEGG IDs for GSEA
#> Created rankings for 243 KEGG metabolites using 'signed_pval' method
#> Ranking range: -10.812 to 8.832
#> Preparing GMT data from 84 pathways
#> GMT preparation summary:
#>   Total pathways processed: 84
#>   Pathways kept (size >= 5): 40
#>   Pathways removed (size < 5): 44
#>   Pathway size distribution:
#>     Min: 5
#>     Median: 7
#>     Max: 21
#>     Mean: 8.9
#> Testing 40 pathways with GSEA
#> GSEA completed: 40 pathways tested, 1 significant at FDR < 0.05
#> Top pathway: 'Aminoacyl-tRNA biosynthesis' (NES = -1.86, FDR = 0.026)
#> Applied KEGG pathway name mapping using 'kegg_id' and 'name' columns
#> Running centrality analysis...
#> Generating metabolite-pathway network visualization...
#> Using top 10 metabolites by centrality
#> Generating enrichment heatmap...
#> Selected top 20 metabolites by enrichment significance for heatmap
#> Generating pathway membership plot...
#> === MEMBERSHIP PLOT DEBUG ===
#> Input metabolites: 94
#> Top_n parameter: 20
#> Initial matched metabolites: 72
#> Filtering to top 20 metabolites by centrality...
#> Total metabolites with centrality: 2961
#> Top metabolites selected: 20
#> Sample top metabolites: C00022, C00026, C00103, C00082, C00249, C00301
#> After centrality filtering - unique metabolites: 20
#> Matrix before occurrence filtering: 20 x 51
#> Pathway counts range: 1 to 4
#> Metabolite counts range: 2 to 24
#> Pathways to keep: 30/51
#> Metabolites to keep: 20/20
#> Final matrix dimensions: 20 x 30
#> === END DEBUG ===
#> Generating Reactome interaction network...
#> Extracting KEGG IDs from input metabolites...
#> Successfully extracted 94 KEGG IDs from 94 input metabolites
#> Sample extracted KEGG IDs: C00209, C00249, C06424, C01727, C04025, C00366
#> Found 25 metabolites with 292 associated Reactome reactions
#> Building metabolite co-occurrence edges from shared reactions...
#> Found 31 valid metabolite pairs sharing >= 1 reaction(s)
#> Applied KEGG pathway name mapping
#> Display names sample: Pyruvate, 2-Oxoglutarate, Uracil, Adenosine, Guanine, Butanoic acid
#> Creating graph with 15 vertices
#> Vertex attributes: name, display_name
#> Graph vertex attributes: name, display_name
#> Sample vertex display_names: Pyruvate, 2-Oxoglutarate, Uracil, Adenosine, Guanine, Butanoic acid
#> Using layout: gem (spacing score: 419.2)

A brief comparison of top pathways:

head(results$pathway_enrichment_results[, c("Pathway", "P_value", "Adjusted_P_value")], 5)
#>                              Pathway      P_value Adjusted_P_value
#> 1              Pyrimidine metabolism 2.729414e-06     0.0001637648
#> 2 Taurine and hypotaurine metabolism 1.598102e-05     0.0004506584
#> 3                  Purine metabolism 2.253292e-05     0.0004506584
#> 4                Caffeine metabolism 1.412711e-04     0.0021190658
#> 5            beta-Alanine metabolism 6.915936e-03     0.0829912279
head(results_bg$pathway_enrichment_results[, c("Pathway", "P_value", "Adjusted_P_value")], 5)
#>                                     Pathway    P_value Adjusted_P_value
#> 1                     Pyrimidine metabolism 0.01369208        0.5172570
#> 2                       Caffeine metabolism 0.01724190        0.5172570
#> 3        Taurine and hypotaurine metabolism 0.04092823        0.8185645
#> 4                         Purine metabolism 0.07233667        0.9612817
#> 5 Biosynthesis of various nucleotide sugars 0.11179352        0.9612817

Network centrality

calculate_metabolite_centrality() ranks metabolites from the pathway–metabolite annotation graph only. Metabolites that share pathways are linked; relative betweenness centrality (RBC) measures how often a metabolite lies on paths between others. Higher RBC means a more bridging role in the pathway map. Those may only decide which metabolites are selected as input; once the set is fixed, ranking is purely topological.

Visualisations

results$pathway_plot

Figure 1. Pathway enrichment. Points correspond to pathways; the horizontal axis shows \(-\log_{10}\) p-value from Fisher’s exact test. Point size reflects the number of overlapping metabolites.

results$impact_plot

Figure 2. Pathway impact. Each pathway is placed by enrichment significance (\(-\log_{10}\) p-value) and a topology-based impact score, so pathways that are both statistically supported and structurally central in the network stand out.

results$gsea_plot

Figure 3. MetSEA (fgsea). Pathways are ranked by enrichment of a metabolite list ordered by statistical evidence and effect size.

results$rbc_plot

Figure 4. Relative betweenness centrality of query metabolites within the pathway–metabolite network (topology only).

results$network_plot

Figure 5. Pathway–metabolite network for top enriched pathways and their overlapping query metabolites.

results$heatmap_plot

Figure 6. Pathway Heatmap. Metabolites are displayed on the y-axis and pathways on the x-axis. The heatmap uses a gradient of −log10 p-values to indicate the significance of enriched pathways and the involvement of each metabolite within the corresponding pathway.

ComplexHeatmap::draw(results$membership_plot)

Figure 7. Pathway Membership plot. Metabolites are displayed on the y-axis and pathways on the x-axis. The heatmap plot illustrates pathway membership by showing each pathway along with its associated metabolites.Only pathways with at least two mapped metabolites are included.

results$interaction_plot

Figure 8. Metabolite interaction network based on shared Reactome reactions (runtime KEGG–ChEBI–Reactome mapping). Edges indicate co-occurrence in the same reaction; node size scales with degree.

Differential analysis followed by enrichment

For an end-to-end analysis that starts from the metabolomics matrix, first run differential analysis with run_de(), then pass the result to enrichmet() via da_results. Significant metabolites are taken from the differential analysis output, so inputMetabolites can be left NULL.

Prepare the abundance matrix in the same format as the example data set metabolomics_mat: metabolites as rows, samples as columns, with clear group labels for the contrasts passed to run_de().

metabolomics_path <- system.file(
    "extdata", "example_data.csv",
    package = "enrichmet"
)

if (metabolomics_path == "") {
    stop("Example file 'example_data.csv' not found in inst/extdata/")
}

metabolomics_mat <- read.csv(
    metabolomics_path,
    row.names = 1,
    check.names = FALSE
)
da_out <- run_de(
    metabolomics_mat,
    "TK-CMV",
    "K-CMV",
    fc_threshold = 1,
    pval_threshold = 0.05
)
#> Volcano plot created successfully
#> Total significant metabolites: 106
#>   Up-regulated: 42
#>   Down-regulated: 64
#> Differential analysis completed successfully
#> Total metabolites analyzed: 345
#> Significant metabolites (FDR < 0.05, |FC| > 1): 106

results_da <- enrichmet(
    inputMetabolites = NULL,
    PathwayVsMetabolites = PathwayVsMetabolites,
    da_results = da_out,
    example_data = example_data,
    kegg_lookup = kegg_lookup,
    reactome_df = reactome_df,
    top_n = 20,
    p_value_cutoff = 1,
    network_top_n = 10,
    heatmap_top_n = 20,
    membership_top_n = 20,
    min_pathway_occurrence = 2,
    min_metabolite_occurrence = 1
)
#> Using metabolites from da_results
#> Using existing 'Significant' column from DA results
#> Using 106 significant metabolites (both up and down) from existing Significant column
#> Processing complex KEGG IDs (splitting by |)...
#> Split 86 complex IDs into 86 unique KEGG IDs
#> Including volcano plot from differential analysis results
#> Running pathway enrichment analysis...
#> Enrichment analysis completed: 59 pathways tested, 59 pathways passed filtering (p <= 1.000000)
#> Running GSEA analysis...
#> Extracting KEGG IDs from 'met_id' column...
#> Successfully extracted KEGG IDs for 398 metabolites
#> Sample extracted KEGG IDs: C00245, C01879, C00025, C00209, C00158, C00049
#> Prepared 243 metabolites with KEGG IDs for GSEA
#> Created rankings for 243 KEGG metabolites using 'signed_pval' method
#> Ranking range: -10.812 to 8.832
#> Preparing GMT data from 84 pathways
#> GMT preparation summary:
#>   Total pathways processed: 84
#>   Pathways kept (size >= 5): 40
#>   Pathways removed (size < 5): 44
#>   Pathway size distribution:
#>     Min: 5
#>     Median: 7
#>     Max: 21
#>     Mean: 8.9
#> Testing 40 pathways with GSEA
#> GSEA completed: 40 pathways tested, 1 significant at FDR < 0.05
#> Top pathway: 'Aminoacyl-tRNA biosynthesis' (NES = -1.85, FDR = 0.036)
#> Applied KEGG pathway name mapping using 'kegg_id' and 'name' columns
#> Running centrality analysis...
#> Generating metabolite-pathway network visualization...
#> Using top 10 metabolites by centrality
#> Generating enrichment heatmap...
#> Selected top 20 metabolites by enrichment significance for heatmap
#> Generating pathway membership plot...
#> === MEMBERSHIP PLOT DEBUG ===
#> Input metabolites: 86
#> Top_n parameter: 20
#> Initial matched metabolites: 66
#> Filtering to top 20 metabolites by centrality...
#> Total metabolites with centrality: 2961
#> Top metabolites selected: 20
#> Sample top metabolites: C00022, C00026, C00082, C00249, C00301, C00219
#> After centrality filtering - unique metabolites: 20
#> Matrix before occurrence filtering: 20 x 51
#> Pathway counts range: 1 to 5
#> Metabolite counts range: 2 to 24
#> Pathways to keep: 27/51
#> Metabolites to keep: 20/20
#> Final matrix dimensions: 20 x 27
#> === END DEBUG ===
#> Generating Reactome interaction network...
#> Extracting KEGG IDs from input metabolites...
#> Successfully extracted 86 KEGG IDs from 86 input metabolites
#> Sample extracted KEGG IDs: C00209, C00249, C06424, C01727, C04025, C00366
#> Found 23 metabolites with 281 associated Reactome reactions
#> Building metabolite co-occurrence edges from shared reactions...
#> Found 29 valid metabolite pairs sharing >= 1 reaction(s)
#> Applied KEGG pathway name mapping
#> Display names sample: Pyruvate, 2-Oxoglutarate, Uracil, Adenosine, Guanine, Uridine
#> Creating graph with 13 vertices
#> Vertex attributes: name, display_name
#> Graph vertex attributes: name, display_name
#> Sample vertex display_names: Pyruvate, 2-Oxoglutarate, Uracil, Adenosine, Guanine, Uridine
#> Using layout: gem (spacing score: 395.28)

When da_results is supplied, metabolites are by default restricted to those that meet a statistical significance threshold (e.g. adjusted p-value). Optionally, users can select metabolites of interest using custom fold-change and p-value cutoffs from the summary statistics, or pass an explicit inputMetabolites vector.

Using a SummarizedExperiment as input

enrichmet() offers an optional, alternative way to supply differential results. Users who already keep their metabolomics data in a SummarizedExperiment do not need to extract a standalone statistics table before running pathway analysis. Instead, the required columns can be read directly from rowData() and passed to the existing da_results interface. This keeps the differential statistics in a single location, the SummarizedExperiment, while enrichmet() continues to use its standard da_results argument for enrichment, MetSEA, centrality, and visualisation. Users who prefer to pass a plain da_out object from run_de() can continue to do so; the two paths are fully interchangeable.

The following rowData() columns are used:

met_id is the metabolite identifier (KEGG ID, name, or mixed; governed by input_type);

log2fc is the log2 fold change, used for MetSEA ranking;

pval or padj provide statistical significance, used to select the query set;

optionally ave_expr and Significant, preserved for downstream use.

library(SummarizedExperiment)
#> Loading required package: MatrixGenerics
#> Loading required package: matrixStats
#> 
#> Attaching package: 'MatrixGenerics'
#> The following objects are masked from 'package:matrixStats':
#> 
#>     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
#>     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
#>     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
#>     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
#>     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
#>     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
#>     colWeightedMeans, colWeightedMedians, colWeightedSds,
#>     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
#>     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
#>     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
#>     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
#>     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
#>     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
#>     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
#>     rowWeightedSds, rowWeightedVars
#> Loading required package: GenomicRanges
#> Loading required package: stats4
#> Loading required package: BiocGenerics
#> Loading required package: generics
#> 
#> Attaching package: 'generics'
#> The following objects are masked from 'package:base':
#> 
#>     as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
#>     setequal, union
#> 
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:stats':
#> 
#>     IQR, mad, sd, var, xtabs
#> The following object is masked from 'package:utils':
#> 
#>     data
#> The following objects are masked from 'package:base':
#> 
#>     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
#>     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
#>     get, grep, grepl, is.unsorted, lapply, Map, mapply, match, mget,
#>     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
#>     rbind, Reduce, rownames, sapply, saveRDS, scale, sequence, table,
#>     tapply, transform, unique, unsplit, which.max, which.min
#> Loading required package: S4Vectors
#> 
#> Attaching package: 'S4Vectors'
#> The following object is masked from 'package:utils':
#> 
#>     findMatches
#> The following objects are masked from 'package:base':
#> 
#>     expand.grid, I, unname
#> Loading required package: IRanges
#> Loading required package: Seqinfo
#> Loading required package: Biobase
#> Welcome to Bioconductor
#> 
#>     Vignettes contain introductory material; view with
#>     'browseVignettes()'. To cite Bioconductor, see
#>     'citation("Biobase")', and for packages 'citation("pkgname")'.
#> 
#> Attaching package: 'Biobase'
#> The following object is masked from 'package:MatrixGenerics':
#> 
#>     rowMedians
#> The following objects are masked from 'package:matrixStats':
#> 
#>     anyMissing, rowMedians

# -------------------------------------------------------------------
# 1. Example abundance matrix
# -------------------------------------------------------------------
metabolomics_path <- system.file(
    "extdata", "example_data.csv",
    package = "enrichmet"
)

metabolomics_mat <- read.csv(
    metabolomics_path,
    row.names = 1,
    check.names = FALSE
)

# -------------------------------------------------------------------
# 2. Differential analysis
# -------------------------------------------------------------------
da_out <- run_de(
    metabolomics_mat,
    "TK-CMV",
    "K-CMV",
    fc_threshold   = 1,
    pval_threshold = 0.05
)
#> Volcano plot created successfully
#> Total significant metabolites: 106
#>   Up-regulated: 42
#>   Down-regulated: 64
#> Differential analysis completed successfully
#> Total metabolites analyzed: 345
#> Significant metabolites (FDR < 0.05, |FC| > 1): 106

# -------------------------------------------------------------------
# 3. Store differential statistics in a SummarizedExperiment
# -------------------------------------------------------------------
kept     <- da_out$full_results$met_id
mat_kept <- metabolomics_mat[kept, , drop = FALSE]

se <- SummarizedExperiment(
    assays  = list(counts = mat_kept),
    rowData = DataFrame(da_out$full_results)
)

# -------------------------------------------------------------------
# 4. Rebuild da_results from the SummarizedExperiment alone
# -------------------------------------------------------------------
da_out_from_se <- da_results_from_se(
    se,
    fc_threshold      = 1,
    pval_threshold    = 0.05,
    split_complex_ids = TRUE
)

str(da_out_from_se)
#> List of 3
#>  $ full_results :'data.frame':   345 obs. of  6 variables:
#>   ..$ met_id     : chr [1:345] "neg_00021_C00245" "neg_00083_C01879" "neg_00096_C00025|C00979" "neg_00123_C00209" ...
#>   ..$ log2fc     : num [1:345] -0.521 -0.608 -0.326 1.562 0.601 ...
#>   ..$ pval       : num [1:345] 1.73e-05 1.33e-02 8.46e-03 2.89e-02 5.36e-04 ...
#>   ..$ padj       : num [1:345] 0.000149 0.024348 0.01621 0.047741 0.001607 ...
#>   ..$ ave_expr   : num [1:345] 26.5 22.5 25 19.5 24 ...
#>   ..$ Significant: chr [1:345] "Not significant" "Not significant" "Not significant" "Up" ...
#>  $ kegg_ready   :'data.frame':   406 obs. of  7 variables:
#>   ..$ met_id     : chr [1:406] "neg_00021_C00245" "neg_00083_C01879" "neg_00096_C00025|C00979" "neg_00096_C00025|C00979" ...
#>   ..$ log2fc     : num [1:406] -0.521 -0.608 -0.326 -0.326 1.562 ...
#>   ..$ pval       : num [1:406] 1.73e-05 1.33e-02 8.46e-03 8.46e-03 2.89e-02 ...
#>   ..$ padj       : num [1:406] 0.000149 0.024348 0.01621 0.01621 0.047741 ...
#>   ..$ ave_expr   : num [1:406] 26.5 22.5 25 25 19.5 ...
#>   ..$ Significant: chr [1:406] "Not significant" "Not significant" "Not significant" "Not significant" ...
#>   ..$ kegg_id    : chr [1:406] "C00245" "C01879" "C00025" "C00979" ...
#>  $ summary_stats:List of 8
#>   ..$ total_metabolites      : int 345
#>   ..$ significant_metabolites: int 106
#>   ..$ upregulated            : int 42
#>   ..$ downregulated          : int 64
#>   ..$ fc_threshold           : num 1
#>   ..$ pval_threshold         : num 0.05
#>   ..$ split_complex_ids      : logi TRUE
#>   ..$ unique_kegg_ids        : int 290

# -------------------------------------------------------------------
# 5. Run enrichmet using only the SE derived object
# -------------------------------------------------------------------
PathwayVsMetabolites <- fetch_kegg_pathway_metabolites(organism = "hsa")
#> Filter applied: excluded 280 pathways. Remaining: 92
kegg_lookup          <- fetch_kegg_compound_lookup()
reactome_df          <- fetch_kegg_reactome(
    kegg_ids = extract_kegg_ids(PathwayVsMetabolites),
    species  = "Homo sapiens"
)
#> Pathway universe contains 2960 unique KEGG compound IDs

results_da <- enrichmet(
    inputMetabolites          = NULL,
    PathwayVsMetabolites      = PathwayVsMetabolites,
    da_results                = da_out_from_se,
    example_data              = example_data,
    kegg_lookup               = kegg_lookup,
    reactome_df               = reactome_df,
    top_n                     = 20,
    p_value_cutoff            = 1,
    network_top_n             = 10,
    heatmap_top_n             = 20,
    membership_top_n          = 20,
    min_pathway_occurrence    = 2,
    min_metabolite_occurrence = 1
)
#> Using metabolites from da_results
#> Using existing 'Significant' column from DA results
#> Using 116 significant metabolites (both up and down) from existing Significant column
#> Processing complex KEGG IDs (splitting by |)...
#> Split 94 complex IDs into 94 unique KEGG IDs
#> Running pathway enrichment analysis...
#> Enrichment analysis completed: 60 pathways tested, 60 pathways passed filtering (p <= 1.000000)
#> Running GSEA analysis...
#> Extracting KEGG IDs from 'met_id' column...
#> Successfully extracted KEGG IDs for 398 metabolites
#> Sample extracted KEGG IDs: C00245, C01879, C00025, C00209, C00158, C00049
#> Prepared 243 metabolites with KEGG IDs for GSEA
#> Created rankings for 243 KEGG metabolites using 'signed_pval' method
#> Ranking range: -10.812 to 8.832
#> Preparing GMT data from 84 pathways
#> GMT preparation summary:
#>   Total pathways processed: 84
#>   Pathways kept (size >= 5): 40
#>   Pathways removed (size < 5): 44
#>   Pathway size distribution:
#>     Min: 5
#>     Median: 7
#>     Max: 21
#>     Mean: 8.9
#> Testing 40 pathways with GSEA
#> GSEA completed: 40 pathways tested, 0 significant at FDR < 0.05
#> Applied KEGG pathway name mapping using 'kegg_id' and 'name' columns
#> Running centrality analysis...
#> Generating metabolite-pathway network visualization...
#> Using top 10 metabolites by centrality
#> Generating enrichment heatmap...
#> Selected top 20 metabolites by enrichment significance for heatmap
#> Generating pathway membership plot...
#> === MEMBERSHIP PLOT DEBUG ===
#> Input metabolites: 94
#> Top_n parameter: 20
#> Initial matched metabolites: 72
#> Filtering to top 20 metabolites by centrality...
#> Total metabolites with centrality: 2961
#> Top metabolites selected: 20
#> Sample top metabolites: C00022, C00026, C00103, C00082, C00249, C00301
#> After centrality filtering - unique metabolites: 20
#> Matrix before occurrence filtering: 20 x 51
#> Pathway counts range: 1 to 4
#> Metabolite counts range: 2 to 24
#> Pathways to keep: 30/51
#> Metabolites to keep: 20/20
#> Final matrix dimensions: 20 x 30
#> === END DEBUG ===
#> Generating Reactome interaction network...
#> Extracting KEGG IDs from input metabolites...
#> Successfully extracted 94 KEGG IDs from 94 input metabolites
#> Sample extracted KEGG IDs: C00209, C00249, C06424, C01727, C04025, C00366
#> Found 25 metabolites with 292 associated Reactome reactions
#> Building metabolite co-occurrence edges from shared reactions...
#> Found 31 valid metabolite pairs sharing >= 1 reaction(s)
#> Applied KEGG pathway name mapping
#> Display names sample: Pyruvate, 2-Oxoglutarate, Uracil, Adenosine, Guanine, Butanoic acid
#> Creating graph with 15 vertices
#> Vertex attributes: name, display_name
#> Graph vertex attributes: name, display_name
#> Sample vertex display_names: Pyruvate, 2-Oxoglutarate, Uracil, Adenosine, Guanine, Butanoic acid
#> Using layout: gem (spacing score: 403.31)

The SummarizedExperiment interface is intended for users who already keep their metabolomics data in a Bioconductor container and want to avoid converting to and from plain data frames. The same input_type options (“kegg”, “name”, “mixed”) apply to the identifiers found in rowData().

Application to a lipid pathway map

The same over-representation interface accepts any pathway–feature table in Pathway / Metabolites format. Lipid ontology content is not shipped with the package. Instead, LION lipid ontology content is obtained at runtime through the BioPortal / LION programming interface when users call the corresponding package helper, such as fetch_lion_lipid_ontology(). Retrieved content is fetched on demand and cached locally (for example, using BiocFileCache) to avoid unnecessary repeated downloads.

Users may also supply their own pathway–feature mapping when appropriate. The example below uses a minimal illustrative table only.

#Creating a toy lipid ontology map

#pathway_lipids <- fetch_lion_lipid_ontology()

#head(pathway_lipids)

# Create a lipids list
#lipids <- c("PA(18:0/16:0)", "PC(O-14:0/18:0)", "PA(22:4(7Z,10Z,13Z,16Z)/18:0)", "Cer(d18:0/24:0)", "SM(d18:0/24:0)", "PS(18:0/16:0)", "TG(14:0/18:3(6Z,9Z,12Z)/22:6(4Z,7Z,10Z,13Z,16Z,19Z))", "DG(16:0/16:0/0:0)[rac]", "TG(16:0/18:3(9Z,12Z,15Z)/22:6(4Z,7Z,10Z,13Z,16Z,19Z))", "PC(O-16:0/20:3(8Z,11Z,14Z))", "PC(10:0/21:0)", "PE(18:0/16:0)", "PA(16:0/18:1(11Z))", "SM(d18:0/22:0)", "TG(17:0/17:1(9Z)/18:0)", "PC(20:0/16:0)", "PC(O-16:0/18:0)", "PS(18:0/22:6(4Z,7Z,10Z,13Z,16Z,19Z))", "TG(18:1(9Z)/18:2(9Z,12Z)/18:2(9Z,12Z))", "SM(d17:1/24:1)", "TG(16:0/16:0/16:0)", "GlcCer(d18:0/24:0)", "Temnoside B", "TG(18:1(9Z)/18:1(9Z)/18:2(9Z,12Z))", "DG(18:1(9Z)/20:3(8Z,11Z,14Z)/0:0)", "PS(16:0/16:0)", "PG(18:1(9Z)/22:4(7Z,10Z,13Z,16Z))", "TG(18:1(9Z)/18:2(9Z,12Z)/18:3(9Z,12Z,15Z))", "SM(d18:2/22:1)", "PC(13:0/22:2(13Z,16Z))", "PE(20:2(5Z,8Z)/18:0)", "TG(12:0/16:0/18:0)", "SM(d18:2/24:1)", "DG(18:3(9Z,12Z,15Z)/20:1(11Z)/0:0)", "TG(16:1(9Z)/18:3(9Z,12Z,15Z)/22:6(4Z,7Z,10Z,13Z,16Z,19Z))", "PG(18:0/16:0)", "TG(18:3(9Z,12Z,15Z)/18:3(9Z,12Z,15Z)/20:4(5Z,8Z,11Z,14Z))", "PE(22:1(11Z)/20:4(5Z,8Z,11Z,14Z))", "SM(d18:0/20:0)", "PG(16:0/16:0)", "PG(18:0/20:4(5Z,8Z,11Z,14Z))", "DG(18:1(9Z)/18:2(9Z,12Z)/0:0)", "Cer(d18:2/20:0)", "PE(O-20:0/20:5(5Z,8Z,11Z,14Z,17Z))", "PI(19:0/22:1(11Z))", "PE(22:1(11Z)/22:4(7Z,10Z,13Z,16Z))", "PC(17:0/22:4(7Z,10Z,13Z,16Z))", "PG(16:0/18:1(9Z))", "PC(O-20:0/22:6(4Z,7Z,10Z,13Z,16Z,19Z))", "TG(16:1(9Z)/18:2(9Z,12Z)/18:2(9Z,12Z))", "SM(d18:1/19:0)", "PE(O-16:0/22:4(7Z,10Z,13Z,16Z))", "TG(18:0/20:4(5Z,8Z,11Z,14Z)/20:5(5Z,8Z,11Z,14Z,17Z))", "PC(22:1(11Z)/22:4(7Z,10Z,13Z,16Z))", "DG(18:1(9Z)/22:4(7Z,10Z,13Z,16Z)/0:0)", "PS(20:1(11Z)/22:4(7Z,10Z,13Z,16Z))", "PS(22:4(7Z,10Z,13Z,16Z)/18:0)", "TG(16:0/16:0/18:0)", "DG(20:1(11Z)/20:4(5Z,8Z,11Z,14Z)/0:0)", "PC(17:0/18:2(9Z,12Z))", "TG(18:0/18:2(9Z,12Z)/18:2(9Z,12Z))", "TG(16:0/16:0/18:3(9Z,12Z,15Z))")
# Perform the enrichment test
#enrichment_results <- perform_enrichment_analysis(lipids, pathway_lipids, p_value_cutoff=0.05)


#create_enrichment_plot(enrichment_results)

Bibliography

Kanehisa, M., Furumichi, M., Sato, Y., Matsuura, Y., & Ishiguro-Watanabe, M. (2025). KEGG: biological systems database as a model of the real world. Nucleic Acids Research, 53, D672–D677. https://www.kegg.jp/
Milacic, M., Beavers, D., Conley, P., Gong, C., Gillespie, M., Griss, J., Haw, R., Jassal, B., Matthews, L., May, B., Petryszak, R., Ragueneau, E., Rothfels, K., Sevilla, C., Shamovsky, V., Stephan, R., Tiwari, K., Varusai, T., Weiser, J., Wright, A., Wu, G., Stein, L., Hermjakob, H., & D’Eustachio, P. (2024). The Reactome Pathway Knowledgebase 2024. Nucleic Acids Research, 52(D1), D672–D678. https://reactome.org/
Molenaar, M. R., Jeucken, A., Wassenaar, T. A., van de Lest, C. H. A., Brouwers, J. F., & Helms, J. B. (2019). LION/web: a web-based ontology enrichment tool for lipidomic data analysis. GigaScience, 8(6), giz061. http://www.lipidontology.com/

Session information

sessionInfo()
#> R version 4.6.2 (2026-10-29)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] SummarizedExperiment_1.43.0 Biobase_2.73.2             
#>  [3] GenomicRanges_1.65.4        Seqinfo_1.3.2              
#>  [5] IRanges_2.47.5              S4Vectors_0.51.10          
#>  [7] BiocGenerics_0.59.12        generics_0.1.4             
#>  [9] MatrixGenerics_1.25.0       matrixStats_1.5.0          
#> [11] enrichmet_0.99.18           BiocStyle_2.41.0           
#> 
#> loaded via a namespace (and not attached):
#>   [1] DBI_1.3.0             gridExtra_2.3.1       httr2_1.3.0          
#>   [4] rlang_1.3.0           magrittr_2.0.5        clue_0.3-68          
#>   [7] GetoptLong_1.1.1      otel_0.2.0            compiler_4.6.2       
#>  [10] RSQLite_3.53.3        reshape2_1.4.5        png_0.1-9            
#>  [13] vctrs_0.7.3           stringr_1.6.0         pkgconfig_2.0.3      
#>  [16] shape_1.4.6.1         crayon_1.5.3          fastmap_1.2.0        
#>  [19] dbplyr_2.6.0          XVector_0.53.0        labeling_0.4.3       
#>  [22] ggraph_2.2.2          rmarkdown_2.32        purrr_1.2.2          
#>  [25] bit_4.6.0             xfun_0.62             cachem_1.1.0         
#>  [28] jsonlite_2.0.0        blob_1.3.0            DelayedArray_0.39.8  
#>  [31] BiocParallel_1.47.0   tweenr_2.0.3          parallel_4.6.2       
#>  [34] cluster_2.1.8.3       R6_2.6.1              bslib_0.12.0         
#>  [37] stringi_1.8.9         RColorBrewer_1.1-3    limma_3.99.4         
#>  [40] jquerylib_0.1.4       Rcpp_1.1.2            iterators_1.0.14     
#>  [43] knitr_1.52            splines_4.6.2         Matrix_1.7-6         
#>  [46] igraph_2.3.4          tidyselect_1.2.1      qvalue_2.45.0        
#>  [49] abind_1.4-8           yaml_2.3.12           viridis_0.6.5        
#>  [52] doParallel_1.0.17     codetools_0.2-20      curl_8.1.0           
#>  [55] plyr_1.8.9            lattice_0.23-1        tibble_3.3.1         
#>  [58] withr_3.0.3           S7_0.2.2              evaluate_1.0.5       
#>  [61] polyclip_1.10-7       BiocFileCache_3.3.1   circlize_0.4.18      
#>  [64] pillar_1.11.1         BiocManager_1.30.27   filelock_1.0.3       
#>  [67] foreach_1.5.2         ggplot2_4.0.3         scales_1.4.0         
#>  [70] glue_1.8.1            maketools_1.3.2       tools_4.6.2          
#>  [73] data.table_1.18.6.1   sys_3.4.3             fgsea_1.39.2         
#>  [76] buildtools_1.0.0      graphlayouts_1.2.5    fastmatch_1.1-8      
#>  [79] cowplot_1.2.0         tidygraph_1.3.1       grid_4.6.2           
#>  [82] tidyr_1.3.2           colorspace_2.1-4      ggforce_0.5.0        
#>  [85] cli_3.6.6             S4Arrays_1.13.2       viridisLite_0.4.3    
#>  [88] ComplexHeatmap_2.29.0 dplyr_1.2.1           gtable_0.3.6         
#>  [91] sass_0.4.10           digest_0.6.39         SparseArray_1.13.4   
#>  [94] ggrepel_0.9.8         rjson_0.2.23          farver_2.1.2         
#>  [97] memoise_2.0.1         htmltools_0.5.9       lifecycle_1.0.5      
#> [100] GlobalOptions_0.1.4   statmod_1.5.2         bit64_4.8.6          
#> [103] MASS_7.3-66

enrichmet: Pathway and Network Analysis for Metabolomics Tutorial