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

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

## ----libraries----------------------------------------------------------------
library(enrichmet)

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

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

## ----plot-pathway, fig.width=10, fig.height=7---------------------------------
results$pathway_plot

## ----plot-impact, fig.width=10, fig.height=7----------------------------------
results$impact_plot

## ----plot-gsea, fig.width=10, fig.height=7------------------------------------
results$gsea_plot

## ----plot-rbc, fig.width=10, fig.height=7-------------------------------------
results$rbc_plot

## ----plot-network, fig.width=10, fig.height=7---------------------------------
results$network_plot

## ----plot-heatmap, fig.width=10, fig.height=7---------------------------------
results$heatmap_plot

## ----plot-membership, fig.width=12, fig.height=8------------------------------
ComplexHeatmap::draw(results$membership_plot)

## ----plot-interaction, fig.width=10, fig.height=8-----------------------------
results$interaction_plot

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

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

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


## ----session-info-------------------------------------------------------------
sessionInfo()