## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

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

## ---- warning = FALSE, message = FALSE, comment = FALSE-----------------------
library(fobitools)

## ---- warning = FALSE, message = FALSE, comment = FALSE-----------------------
# CRAN
library(tidyverse)
library(ggrepel)
library(kableExtra)

# Bioconductor
library(POMA)
library(metabolomicsWorkbenchR)
library(SummarizedExperiment)

## ---- warning = FALSE, message = FALSE, comment = FALSE-----------------------
data <- do_query(
  context = "study",
  input_item = "analysis_id",
  input_value = "AN000961",
  output_item = "SummarizedExperiment")

## ---- warning = FALSE, message = FALSE, comment = FALSE-----------------------
## features 
features <- assay(data)
rownames(features) <- rowData(data)$metabolite

## metadata
pdata <- colData(data) %>%
  as.data.frame() %>%
  tibble::rownames_to_column("ID") %>%
  relocate(grouping, .before = Sodium.level) %>%
  mutate(grouping = case_when(grouping == "A(salt sensitive)" ~ "A",
                              grouping == "B(salt insensitive)" ~ "B"))

## ---- warning = FALSE, message = FALSE, comment = FALSE-----------------------
data_sumexp <- PomaSummarizedExperiment(target = pdata, features = t(features))

## ---- warning = FALSE, message = FALSE, comment = FALSE-----------------------
data_preprocessed <- data_sumexp %>%
  PomaImpute(ZerosAsNA = TRUE, cutoff = 20, method = "knn") %>%
  PomaNorm(method = "log_pareto") %>%
  PomaOutliers(coef = 3)

## ---- warning = FALSE, message = FALSE, comment = FALSE-----------------------
limma_res <- data_preprocessed %>%
  PomaLimma(contrast = "A-B", adjust = "fdr") %>%
  dplyr::rename(ID = feature)

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

## ---- warning = FALSE, message = FALSE, comment = FALSE, eval = FALSE---------
#  cat(limma_res$ID, sep = "\n")

## ----MAconvertID, message = FALSE, warning = FALSE, comment = FALSE, echo = FALSE, fig.align = "center", out.width = "100%", fig.cap = 'Metabolite names conversion using MetaboAnalyst.'----
knitr::include_graphics("figure/MetaboAnalyst_convertID.png")

## ---- warning = FALSE, message = FALSE, comment = FALSE, eval = FALSE---------
#  ST000629_MetaboAnalyst_MapNames <- readr::read_delim("https://www.metaboanalyst.ca/MetaboAnalyst/resources/users/XXXXXXX/name_map.csv", delim = ",")

## ---- warning = FALSE, message = FALSE, comment = FALSE, echo = FALSE---------
load("data/ST000629_MetaboAnalyst_MapNames.rda")

## ---- warning = FALSE, message = FALSE, comment = FALSE-----------------------
annotated_limma <- ST000629_MetaboAnalyst_MapNames %>%
  dplyr::rename(ID = Query) %>%
  right_join(limma_res, by = "ID")

limma_FOBI_names <- annotated_limma %>%
  dplyr::pull("KEGG") %>%
  fobitools::id_convert(to = "FOBI")

# show the ID conversion results
limma_FOBI_names %>%
  head() %>%
  kbl(row.names = FALSE, booktabs = TRUE) %>%
  kable_styling(latex_options = c("striped"))

## ---- warning = FALSE, message = FALSE, comment = FALSE-----------------------
limma_FOBI_names <- limma_FOBI_names %>%
  right_join(annotated_limma, by = "KEGG") %>%
  dplyr::select(FOBI, KEGG, ID, logFC, P.Value, adj.P.Val) %>%
  dplyr::arrange(-dplyr::desc(P.Value))

## ---- warning = FALSE, message = FALSE, comment = FALSE-----------------------
metaboliteList <- limma_FOBI_names$FOBI[limma_FOBI_names$P.Value < 0.05]
metaboliteUniverse <- limma_FOBI_names$FOBI

fobitools::ora(metaboliteList = metaboliteList,
               metaboliteUniverse = metaboliteUniverse,
               pvalCutoff = 1) %>%
  kbl(row.names = FALSE, booktabs = TRUE) %>%
  kable_styling(latex_options = c("striped"))

## ---- warning = FALSE, message = FALSE, comment = FALSE-----------------------
limma_FOBI_msea <- limma_FOBI_names %>%
  select(FOBI, P.Value) %>%
  filter(!is.na(FOBI)) %>%
  dplyr::arrange(-dplyr::desc(abs(P.Value)))

FOBI_msea <- as.vector(limma_FOBI_msea$P.Value)
names(FOBI_msea) <- limma_FOBI_msea$FOBI

msea_res <- fobitools::msea(FOBI_msea, pvalCutoff = 1)

msea_res %>%
  kbl(row.names = FALSE, booktabs = TRUE) %>%
  kable_styling(latex_options = c("striped"))

## ---- warning = FALSE, message = FALSE, comment = FALSE, fig.width = 12, fig.height = 8----
ggplot(msea_res, aes(x = -log10(pval), y = NES, color = NES, size = classSize, label = className)) +
  xlab("-log10(P-value)") +
  ylab("NES (Normalized Enrichment Score)") +
  geom_point() +
  ggrepel::geom_label_repel(color = "black", size = 7) +
  theme_bw() +
  theme(legend.position = "top",
        text = element_text(size = 22)) +
  scale_color_viridis_c() +
  scale_size(guide = "none")

## ---- warning = FALSE, message = FALSE, comment = FALSE, fig.width = 12, fig.height = 10----
FOBI_terms <- msea_res %>%
  unnest(cols = leadingEdge)

fobitools::fobi %>%
  filter(FOBI %in% FOBI_terms$leadingEdge) %>%
  pull(id_code) %>%
  fobi_graph(get = "anc",
             labels = TRUE,
             legend = TRUE,
             labelsize = 6,
             legendSize = 20)

## -----------------------------------------------------------------------------
sessionInfo()