---
title: "enrichmet: Pathway and Network Analysis for Metabolomics Tutorial"
author: "Yonatan Ayalew Mekonnen"
date: "`r Sys.Date()`"
output:
  BiocStyle::html_document:
    toc: true
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{enrichmet: Pathway and Network Analysis for Metabolomics}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    fig.align = "center",
    message = FALSE,
    warning = FALSE
)
```

# 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

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

```{r libraries}
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:

1. `inputMetabolites` — character vector of KEGG IDs and/or names (see
   Section 3.3); or
2. `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):

```{r pathway-map, eval=FALSE}
# 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.

```{r case-study, message=TRUE}
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")
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)
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"
)
```

## 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.

```{r background, message=TRUE}
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"
)
```

A brief comparison of top pathways:

```{r compare-bg}
head(results$pathway_enrichment_results[, c("Pathway", "P_value", "Adjusted_P_value")], 5)
head(results_bg$pathway_enrichment_results[, c("Pathway", "P_value", "Adjusted_P_value")], 5)
```

# 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

```{r plot-pathway, fig.width=10, fig.height=7}
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.

```{r plot-impact, fig.width=10, fig.height=7}
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.

```{r plot-gsea, fig.width=10, fig.height=7}
results$gsea_plot
```

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

```{r plot-rbc, fig.width=10, fig.height=7}
results$rbc_plot
```

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

```{r plot-network, fig.width=10, fig.height=7}
results$network_plot
```

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

```{r plot-heatmap, fig.width=10, fig.height=7}
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.

```{r plot-membership, fig.width=12, fig.height=8}
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.

```{r plot-interaction, fig.width=10, fig.height=8}
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()`.

```{r de-workflow, message=TRUE}
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
)

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
)
```

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.

```{r summarizedExperiment, message=TRUE}
library(SummarizedExperiment)

# -------------------------------------------------------------------
# 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
)

# -------------------------------------------------------------------
# 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)

# -------------------------------------------------------------------
# 5. Run enrichmet using only the SE derived object
# -------------------------------------------------------------------
PathwayVsMetabolites <- fetch_kegg_pathway_metabolites(organism = "hsa")
kegg_lookup          <- fetch_kegg_compound_lookup()
reactome_df          <- fetch_kegg_reactome(
    kegg_ids = extract_kegg_ids(PathwayVsMetabolites),
    species  = "Homo sapiens"
)

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
)
```
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.


```{r lipidomics, message=TRUE}
#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

1. 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/

2. 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/

3. 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

```{r session-info}
sessionInfo()
```