```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r}
suppressPackageStartupMessages({

library(readxl)
library(readr)

library(limma)
library(stringr)
library(rvest)
library(ggrepel)
library(kableExtra)
library(tidyverse)


library(POMA)
library(metabolomicsWorkbenchR)
library(SummarizedExperiment)

library(BiocFileCache)
library(enrichmet)
})

```


# Introduction

This document contains the preprocessing and analysis of two metabolomics datasets that will be used for benchmarking several enrichment analysis packages.

The datasets are:

- The Cachexia dataset a well known eMPLE DATASET USED IN mETABOaNALYST
- The MST000291 dataset used by the fobitools package as an example and available in MetabolomicsWorkbench

# The cachexia dataset

This is an example introduced by Xuia et alt. in MetaboAnalyst and is used in most tutorials of this software.

The dataset is available as a .csv file in the metaboData repository in github (https://github.com/nutrimetabolomics/metaboData/blob/main/Datasets/2024-Cachexia/human_cachexia.csv). For implicity it has been downladed from there and re-stored as an Excel File.

```{r}
human_cachexia <- read_excel("data/human_cachexia.xlsx")
str(human_cachexia)
```

## Selection of differentially quantified metabolites

This dataset contains only metabolites that have been shown to be significantly different among the groups.

This is confirmed by the analysis performed in the following.

```{r}

# 1. Definir el vector de grups a partir dels noms de les columnes
sample_names <- colnames(human_cachexia)[-1]   # excepte 'Metabolite'

group <- ifelse(str_detect(sample_names, "^ctrl"), "Control", "Cachexia")
group <- factor(group)

# 2. Construir la matriu d'expressió (metabolits x mostres)
X <- as.matrix(human_cachexia[ , -1])
rownames(X) <- human_cachexia$Metabolite

# 3. Disseny experimental
design <- model.matrix(~ 0+group)

# 4. Ajustar model amb limma
fit <- lmFit(X, design)
fit <- eBayes(fit)

# 5. Resultats
res <- topTable(fit, n=Inf)
tail(res)

# 6. Mostrar els primers resultats
print(res)
topCachexia <- res; rm(res)
save(topCachexia, file="data/Cachexia_topTable.Rda")
write.table(rownames(res), file="data/cachexia_ids.txt", row.names = FALSE, quote=FALSE)
```

## Selected Metabolites and its identifiers

Metabolites in the table are identified by their common names, a type of identifier very useful for scientists working with these metabolites at chemical or nutritional level. However, most enrichment analysis tools require a more objective ID from public databases such as PubChem, HumanDB or KEGG. The associated identifiers for these metabolites, when available, have been obtained using [MetaboAnalyst Compound ID converter](https://www.metaboanalyst.ca/MetaboAnalyst/upload/ConvertView.xhtml).

```{r}
cachexia_map <- read.csv("data/Cachexia_map.csv")
head(cachexia_map)
```

The obtained identifiers may be used as input for enrichment analysis.

# The ST000291 dataset

This dataset comes from a study that has bee uploaded to Metabolomics Workbench with ID ST0000291: [LC-MS Based Approaches to Investigate Metabolomic Differences in the Urine of Young Women after Drinking Cranberry Juice or Apple Juice](https://www.metabolomicsworkbench.org/data/DRCCMetadata.php?Mode=Study&StudyID=ST000291)

Once downloaded the data consists of three separate files

- The features data
  - 1541 variables, 45 samples

- The metadata
  - 45 rows, two columns (sample name, group name)
  
The metabolite names
  - 1541 row, 3 columns (original name, PubChem ID and KEGG ID)

Indeed the downloaded data are stored in a particular format "MWTab" so some functions to extract data from it have been prepared.

**WARNING**

> Although the original study followed a 2×2 crossover design, the Metabolomics Workbench dataset (ST000291) does not include subject identifiers, and therefore paired or repeated-measures models cannot be applied. As a result, treatment groups must be analysed as independent samples.

> This choice is conservative: ignoring the within-subject pairing increases residual variance and therefore raises the probability of false negatives, not false positives. Consequently, any metabolites that remain significant under this independent-samples analysis are likely to represent strong and biologically robust treatment effects.

## Data loading and preprocessing

Data in Metabolomics Workbench is not preprocessed so some work has to be performed before they can be analyzed.

The preprocessing steps are performed using the POMA Bioconductor package.

```{r}

data_negative_mode <- do_query(
  context = "study",
  input_item = "analysis_id",
  input_value = "AN000465",
  output_item = "SummarizedExperiment")

data_positive_mode <- do_query(
  context = "study",
  input_item = "analysis_id",
  input_value = "AN000464",
  output_item = "SummarizedExperiment")

## ----metabolitenames, message = FALSE, warning = FALSE, comment = FALSE, echo = FALSE, fig.align = "center", out.width = "100%", fig.cap = 'Metabolite identifiers of the ST000291 Metabolomics Workbench study.'----
# knitr::include_graphics("figure/ST000291_names.png")

## ----warning = FALSE, message = FALSE, comment = FALSE------------------------
metaboliteNamesURL <- "https://www.metabolomicsworkbench.org/data/show_metabolites_by_study.php?STUDY_ID=ST000291&SEARCH_TYPE=KNOWN&STUDY_TYPE=MS&RESULT_TYPE=1"

metaboliteNames <- metaboliteNamesURL %>% 
  read_html() %>% 
  html_nodes(".datatable")

metaboliteNames_negative <- metaboliteNames %>%
  .[[1]] %>%
  html_table() %>%
  dplyr::select(`Metabolite Name`, PubChemCompound_ID, `Kegg Id`)

metaboliteNames_positive <- metaboliteNames %>%
  .[[2]] %>%
  html_table() %>%
  dplyr::select(`Metabolite Name`, PubChemCompound_ID, `Kegg Id`)

metaboliteNames <- bind_rows(metaboliteNames_negative, metaboliteNames_positive) %>%
  dplyr::rename(names = 1, PubChem = 2, KEGG = 3) %>%
  mutate(KEGG = ifelse(KEGG == "-", "UNKNOWN", KEGG),
         PubChem = ifelse(PubChem == "-", "UNKNOWN", PubChem)) %>%
  filter(!duplicated(PubChem))

## ----warning = FALSE, message = FALSE, comment = FALSE------------------------
## negative mode features 

features_negative <- assay(data_negative_mode) %>%
  dplyr::slice(-n())
rownames(features_negative) <- rowData(data_negative_mode)$metabolite[1:(length(rowData(data_negative_mode)$metabolite)-1)]

## positive mode features
features_positive <- assay(data_positive_mode) %>%
  dplyr::slice(-n())
rownames(features_positive) <- rowData(data_positive_mode)$metabolite[1:(length(rowData(data_positive_mode)$metabolite)-1)]

## combine positive and negative mode and set PubChem IDs as feature names
features <- bind_rows(features_negative, features_positive) %>%
  tibble::rownames_to_column("names") %>%
  right_join(metaboliteNames, by = "names") %>%
  select(-names, -KEGG) %>%
  tibble::column_to_rownames("PubChem")

## metadata
pdata <- colData(data_negative_mode) %>% # or "data_positive_mode". They are equal
  as.data.frame() %>%
  tibble::rownames_to_column("ID") %>%
  mutate(Treatment = case_when(Treatment == "Baseline urine" ~ "Baseline",
                               Treatment == "Urine after drinking apple juice" ~ "Apple",
                               Treatment == "Urine after drinking cranberry juice" ~ "Cranberry"))

## ----warning = FALSE, message = FALSE, comment = FALSE------------------------
pdata$Treatment <- as.character(pdata$Treatment)
pdata$Treatment <- factor(pdata$Treatment,
                          levels = c("Baseline", "Apple", "Cranberry"))

data_sumexp <- PomaCreateObject(metadata = pdata, features = t(features))


data_preprocessed <- data_sumexp %>%
  PomaImpute() %>%             # sense group_by
  PomaNorm()

```

The normalized, imputed feature matrix is now:

```{r extract-preprocessed}
expr_all <- assay(data_preprocessed)              # rows: features (PubChem), cols: samples
sample_metadata <- colData(data_preprocessed) %>%
  as.data.frame()

dim(expr_all)
head(expr_all[, 1:5])

```

## Selection of differential metabolites

We focus on **Baseline vs Cranberry**, ignoring Apple for the main comparison.

```{r}
limma_resCB <- data_preprocessed %>%
  PomaLimma(contrast = "Cranberry-Baseline", adjust = "fdr") %>%
  dplyr::rename(PubChemCID = feature) %>% 
  dplyr::mutate(PubChemCID = gsub("X", "", PubChemCID)) %>%
  filter(PubChemCID != "UNKNOWN")

# show the first 10 features
limma_resCB %>%
  dplyr::slice(1L:10L) %>%
  kbl(row.names = FALSE, booktabs = TRUE) %>%
  kable_styling(latex_options = c("striped"))

# save(limma_res, file="data/MW_ST000291-limma_res.Rda")
```

## Restictring the analysis to metabolites with HMDB_ID

In order to verify some calculations we wish to compare the results with those obtained with MetaboAnalyst. This requires metabolites to have either KEGG or HMDB identifiers.

Pubchem IDs,converted to other systems are available from file `MW_ST000291__map.csv`

```{r}
MW_map <- readr::read_csv("data/MW_ST000291__map.csv", show_col_types = FALSE) %>%
  janitor::clean_names() %>%
  mutate(
    PubChemCID = as.character(query),  # alinyem nom
    HMDB = ifelse(hmdb == "" | is.na(hmdb), NA, hmdb)
  ) %>%
  select(PubChemCID, HMDB)

head(MW_map)

``` 

First we filter out all metabolites not mapped to HMDB

```{r}
expr_all <- assay(data_preprocessed)              # rows = PubChem, cols = samples
metadata <- colData(data_preprocessed) %>% as.data.frame()

pubchem_ids <- gsub("^X", "", rownames(expr_all))
map_expr <- tibble(
  PubChemCID = pubchem_ids
) %>%
  left_join(MW_map, by = "PubChemCID")

keep <- !is.na(map_expr$HMDB)
data_preprocessed_HMDB <- data_preprocessed[keep, ]
head(assay(data_preprocessed_HMDB))

hmdb_ids <- map_expr$HMDB[keep]
rownames(data_preprocessed)[keep] <- hmdb_ids
data_preprocessed_HMDB <- data_preprocessed[keep, ]
head(assay(data_preprocessed_HMDB))
```

Now the analysis can be reproduced

```{r}
limma_HMDB_CB <- data_preprocessed_HMDB %>%
  PomaLimma(contrast = "Cranberry-Baseline", adjust = "fdr")
head(limma_HMDB_CB)

background_HMDB <- limma_HMDB_CB$feature

selected_HMDB <- limma_HMDB_CB %>%
  filter(adj_pvalue < 0.25) %>%
  pull(feature)

```

The topTables are saved to be used later: 

```{r}
save(limma_resCB, limma_HMDB_CB, file="data/MW_ST000291-limma_res.Rda")
```

We remove unnecessary packages used onlñy in this section

```{r}
detach(package:POMA)
detach(package:SummarizedExperiment)
```


# `enrichmet` package exmple data


The `enrichmet`package has deposited its example data in Zotero. For reproducibility we download the data directly from there:


To run the **enrichmet** package, you need to provide at least the following three files:
  
  - **`PathwayVsMetabolites`**  
  This file defines the mapping between metabolic pathways and their associated metabolites, serving as the background reference for Fisher’s exact test–based enrichment analysis (see section 2.4).

The mappings are derived from the KEGG database (https://www.kegg.jp/), accessed programmatically via the KEGGREST package. To ensure reproducibility and minimize run-time delays, curated human-specific pathway–metabolite mappings are periodically updated and archived on Zenodo [1].
We provide a human-specific `PathwayVsMetabolites` file available here:  
  [Human-PathwaysVsMetabolites.csv](https://zenodo.org/api/records/15498097/files/Human-PathwaysVsMetabolites.csv)  
However, users can supply a file tailored to their species of interest.

- **`InputMetabolites`**  
  A list of metabolites from the user’s experimental dataset for analysis, typically consisting of metabolite identifiers (e.g., significant metabolites from differential analysis (see section 2.4)). This input is not provided by the package and must be prepared by the user. The file serves as the query set for enrichment analysis, and its identifiers should match those used in the background pathway–metabolite mappings (e.g., KEGG compound IDs).

For demonstration purposes, we also include summary statistics derived from an example experimental dataset. These example data are intended solely to showcase the package workflow and will be described in detail in a subsequent publication. 


- **`SummaryStatistics`**  
  The enrichment analysis requires statistical results associated with the input metabolites, such as fold changes, p-values, or other relevant effect size metrics (see section 2.4). These values guide the enrichment procedure by reflecting both the significance and the magnitude of change for each metabolite.

This information is not generated by the package and must be provided by the user. Summary statistics may be supplied either as a standalone data frame or extracted directly from the rowData() of a SummarizedExperiment object, provided that appropriate columns are available for metabolite identifiers, p-values, and fold changes (see section 2.4 example_data).

For demonstration purposes, the package includes an example summary statistics dataset. These example data are provided solely to illustrate the workflow and will be described in detail in a forthcoming publication.

The **enrichmet**() function integrates both p-values and fold changes when performing MetSEA. By considering the direction and magnitude of changes, this approach yields richer biological insights than methods based only on binary classifications of "significant" versus "not significant."

## Optional Input Files

Including these optional files can enhance the analysis and its biological interpretation:
  
  - **`KEGGLookup`**  
  A reference file linking KEGG compound IDs (https://www.kegg.jp/) to metabolite names or other annotations, which helps standardize identifiers and improve pathway mapping [1].
We provide a KEGGLookup file here:  
  [kegg_lookup.xlsx](https://zenodo.org/api/records/15498097/files/kegg_lookup.xlsx)  
Users may also provide a lookup file specific to their species.

- **`STITCHInteractions`**  
  A dataset of known or predicted interactions between metabolites or proteins, obtained from the freely available STITCH database (http://stitch.embl.de/
  ). These interactions are useful for constructing networks to support pathway analysis [2].
We supply an example file here:  
  [chemical_chemical.tsv](https://zenodo.org/api/records/15498097/files/chemical_chemical.tsv)  
Users may provide their own interaction data appropriate for their species.

- **`mapping_df`**  
  A mapping data frame that links STITCH identifiers with metabolite names, facilitating integration between datasets [2].
The provided example file is available here:  
  [mapping_df.xlsx](https://zenodo.org/api/records/15498097/files/mapping_df.xlsx)  
Users can substitute this with a mapping file relevant to their species of interest.


## Case study of KrasG12D mutation

We used summary statistics from differential analyses based on metabolomic measurements comparing KrasG12D mutations with and without TAp73 deletion. The top 50 metabolites, ranked by significance, were selected for pathway enrichment analysis using Fisher’s Exact Test. These summary statistics were also utilized for MetSEA analysis.

```{r}
# InputMetabolites: A vector or list containing the metabolites of interest.
inputMetabolites <- c(
  "C02721", "C03736", "C03139", "C00074", "C02862", "C00092", "C00275",
  "C00386", "C10172", "C00361", "C00242", "C19806", "C16358", "C16357",
  "C16353", "C02140", "C05488", "C01181", "C03017", "C02632", "C00246",
  "C00079", "C08317", "C00588", "C16527", "C00065", "C00249", "C00836",
  "C02990", "C03690", "C04025", "C06424", "C00123", "C00606", "C03067",
  "C00539", "C00633", "C00519", "C00366", "C00148", "C00219", "C08261",
  "C00385", "C06429", "C05637", "C00295", "C06428", "C16513", "C01909",
  "C03299", "C01530", "C01571", "C08278", "C00599", "C15587", "C00144",
  "C03739"
)
```

```{r}
# Standard cache location
cache_dir <- tools::R_user_dir("enrichmet", "cache")
bfc <- BiocFileCache(cache_dir, ask = FALSE)

# Helper function to download & cache any file from Zenodo
get_cached_file <- function(url) {
  bfcrpath(bfc, url)
}

# -------------------------
# Load files with caching
# -------------------------

# PathwayVsMetabolites
pathway_url <- "https://zenodo.org/api/records/15498097/files/Human-PathwaysVsMetabolites.csv/content"
pathway_path <- get_cached_file(pathway_url)
PathwayVsMetabolites <- read.csv(pathway_path)
str(PathwayVsMetabolites)

# SummaryStatistics
example_url <- "https://zenodo.org/api/records/15498097/files/example_data.xlsx/content"
example_path <- get_cached_file(example_url)
example_data <- read_excel(example_path)
show(example_data[1:3,])


# KEGGLookup
kegg_url <- "https://zenodo.org/api/records/15498097/files/kegg_lookup.xlsx/content"
kegg_path <- get_cached_file(kegg_url)
kegg_lookup <- read_excel(kegg_path)
head(kegg_lookup)

# Mapping_df
mapping_url <- "https://zenodo.org/api/records/15498097/files/mapping_df.xlsx/content"
mapping_path <- get_cached_file(mapping_url)
mapping_df <- read_excel(mapping_path)
head(mapping_df)

# STITCHInteractions
# stitch_url <- "https://zenodo.org/records/15498097/files/chemical_chemical.tsv"
# stitch_path <- get_cached_file(stitch_url)
# stitch_df <- read_tsv(stitch_path)
# head(stitch_df)

```


The data from "EnrichMetExample" is not an EnrichmentSet, so a new one is created from the `PathwayvsMetabolites` data.frame:

First we prepare the `long` version of `PathwayVsMetabolites`:

```{r}
tmp <- PathwayVsMetabolites %>%
  mutate(Metabolites = strsplit(Metabolites, ","))

unnested <- tidyr::unnest(tmp, cols = c(Metabolites))

PathwayVsMetabolites_long <- dplyr::rename(unnested, kegg_id = Metabolites)

head(PathwayVsMetabolites_long)

PathwayVsMetabolites_long <- PathwayVsMetabolites_long %>%
  mutate(PathwayID = paste0("PW_", as.integer(as.factor(Pathway))))
head(PathwayVsMetabolites_long)

PathwayVsMetabolites_long <- PathwayVsMetabolites_long %>%
  left_join(kegg_lookup, by = c("kegg_id" = "kegg_id"))
head(PathwayVsMetabolites_long)
```

With this we can directly build the `EnrichmentSet` object:

```{r}
require(localEnrichment)

KEGGset_KEGG <- buildEnrichmentSet(
  data = PathwayVsMetabolites_long,
  id_col = "kegg_id",
  category_col = "Pathway",
  set_id_col = "PathwayID",
  set_name = "KEGG_hsa_pathways",
  source = "KEGG",
  species = "Homo sapiens",
  version = Sys.Date(),
  description = "KEGG human pathways mapped to KEGG compound IDs",
  sep = ";"
)

summary(KEGGset_KEGG)
head(KEGGset_KEGG@data)
```

All the data is saved as a binary file to be used later

```{r}
EM_inputMetabolites <- inputMetabolites
EM_summaryStats <- example_data
EM_PathwayVsMetabolites<- PathwayVsMetabolites
EM_mapping_df <- mapping_df
# EM_stitch_df <- stitch_df
EM_kegg_lookup <- kegg_lookup
save(EM_inputMetabolites, KEGGset_KEGG, EM_PathwayVsMetabolites, EM_summaryStats, EM_kegg_lookup, EM_mapping_df, file="data/enrichmet_example.Rda")

```


