## ----style, echo = FALSE, results = 'asis'------------------------------------
BiocStyle::markdown()

## ----eval=FALSE---------------------------------------------------------------
# if (!require("BiocManager", quietly = TRUE)) {
#   install.packages("BiocManager")
# }
# 
# # The following initializes usage of Bioc devel
# BiocManager::install(version = "devel")
# 
# BiocManager::install("MetaboDynamics")

## ----setup--------------------------------------------------------------------
library(MetaboDynamics)
library(SummarizedExperiment)
library(ggplot2)
library(dplyr)
library(tidyr)

## ----fig.wide=TRUE------------------------------------------------------------
data("longitudinalMetabolomics")
# convert to dataframe
concentrations <- as.data.frame(cbind(
  as.data.frame(t(assays(longitudinalMetabolomics)$concentrations)),
  as.data.frame(colData(longitudinalMetabolomics))
))
concentrations <- concentrations %>% pivot_longer(
  cols = seq_len(4), names_to = "time",
  values_to = "concentration"
)
ggplot(concentrations, aes(x = concentration)) +
  geom_density() +
  theme_bw() +
  facet_grid(cols = vars(time), rows = vars(condition)) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  ggtitle("raw data", "raw measurements")

## ----fig.wide=TRUE------------------------------------------------------------
ggplot(concentrations, aes(x = log(concentration))) +
  geom_density() +
  theme_bw() +
  facet_grid(cols = vars(time), rows = vars(condition)) +
  ggtitle("data", "log-transformed values")

## ----fig.wide=TRUE------------------------------------------------------------
ggplot(concentrations) +
  geom_line(aes(
    x = time, y = log(concentration), col = metabolite,
    group = interaction(metabolite, replicate)
  )) +
  theme_bw() +
  xlab("timepoint") +
  theme(legend.position = "none") +
  facet_grid(rows = vars(condition)) +
  ggtitle("raw metabolite dynamics", "color=metabolite")

## ----fig.wide=TRUE------------------------------------------------------------
data("longitudinalMetabolomics")
# convert to dataframe
concentrations <- as.data.frame(cbind(
  as.data.frame(t(assays(longitudinalMetabolomics)$scaled_log)),
  as.data.frame(colData(longitudinalMetabolomics))
))
concentrations <- concentrations %>% pivot_longer(
  cols = seq_len(4), names_to = "time",
  values_to = "concentration"
)
ggplot(concentrations) +
  geom_line(aes(
    x = time,
    y = concentration, col = metabolite,
    group = interaction(metabolite, replicate)
  )) +
  theme_bw() +
  xlab("timepoint") +
  theme(legend.position = "none") +
  facet_grid(rows = vars(condition)) +
  ggtitle("standardized dynamics", "color=metabolite")

## ----fig.wide=TRUE------------------------------------------------------------
# we can hand a SummarizedExperiment object to the function
data("longitudinalMetabolomics")
# we only use a subsection of the simulated data (1 condition and subsample of
# the whole dataset) for demonstration purposes
samples <- c(
  "UMP", "Taurine", "Succinate", "Phosphocreatine", "PEP",
  "Malic acid", "L-Cystine", "CMP", "Citramalic acid", "2-Aminomuconate"
)
# only one condition
data <- longitudinalMetabolomics[, longitudinalMetabolomics$condition == "A" &
  longitudinalMetabolomics$metabolite %in% samples]

# fit model
data <- fit_dynamics_model(
  data = data, scaled_measurement = "m_scaled",
  assay = "scaled_log", time = "time",
  condition = "condition", max_treedepth = 10,
  adapt_delta = 0.95, # default 0.95
  iter = 5000,
  cores = 1,
  chains = 2 # only set to 2 for vignette, default = 4
)

## -----------------------------------------------------------------------------
# extract diagnostics
data <- diagnostics_dynamics(
  data = data,
  iter = 5000, # number of iterations used for model fitting
  # the dynamic model
  chains = 2 # number of chains used for model fitting
)

plot_diagnostics(
  data = data
)

# PPCs can be accessed with
plot_PPC(
  data = data
)

## ----fig.wide=TRUE------------------------------------------------------------
# #extract estimates
data <- estimates_dynamics(
  data = data, iter = 5000,
  chains = 2, condition = "condition"
)

## ----fig.wide=TRUE------------------------------------------------------------
# 1) the differences between two timepoints
plot_estimates(
  data = data,
  dynamics = FALSE
)

## ----fig.wide=TRUE------------------------------------------------------------
# 2) dynamic profiles
plot_estimates(
  data = data,
  delta_t = FALSE
)

## ----fig.wide=TRUE------------------------------------------------------------
# get distances between vectors
clust_A <- metadata(data)[["estimates_dynamics"]][["A"]]
clust_A <- clust_A %>%
  select(metabolite.ID, condition, time.ID, mu_mean) %>%
  pivot_wider(id_cols = c(metabolite.ID, condition), names_from = time.ID, values_from = mu_mean)

dist <- clust_A %>% select(-c(metabolite.ID, condition))
dd_A <- dist(dist,
  method = "euclidean"
)
# hierarchical clustering
clust <- hclust(dd_A, method = "ward.D2")
clust_cut <- cutree(clust, k = 8)
# adding cluster ID to estimates
clust_A$cluster <- clust_cut
clust_A

## ----fig.wide=TRUE------------------------------------------------------------
data <- metadata(longitudinalMetabolomics)[["cluster"]]
temp <- data
temp <- temp %>% pivot_longer(
  cols = 4:7,
  names_to = "timepoint", values_to = "mu_mean"
)
ggplot(temp, aes(
  x = as.factor(as.numeric(as.factor(timepoint))),
  y = mu_mean, group = metabolite
)) +
  geom_line() +
  xlab("timepoint") +
  ylab("estimated mean concentration") +
  theme_bw() +
  theme(legend.position = "none") +
  facet_grid(rows = vars(condition), cols = vars(cluster)) +
  ggtitle("clustered dynamics", "panels=cluster ID")

## -----------------------------------------------------------------------------
data("metabolite_modules")
help("metabolite_modules")
head(metabolite_modules)
data("modules_compounds")
help("modules_compounds")
head(modules_compounds)

## ----fig.wide=TRUE------------------------------------------------------------
data <- ORA_hypergeometric(
  background = modules_compounds,
  annotations = metabolite_modules,
  data = longitudinalMetabolomics, tested_column = "middle_hierarchy"
)
plot_ORA(data, tested_column = "middle_hierarchy")

data <- ORA_hypergeometric(
  background = modules_compounds,
  annotations = metabolite_modules,
  data = longitudinalMetabolomics,
  tested_column = "lower_hierarchy"
)
plot_ORA(data, tested_column = "lower_hierarchy")

## ----fig.aling='center',fig.dpi=150-------------------------------------------
data("longitudinalMetabolomics")
longitudinalMetabolomics <- compare_dynamics(
  data = longitudinalMetabolomics,
  dynamics = c("mu1_mean", "mu2_mean", "mu3_mean", "mu4_mean"),
  cores = 1
)
heatmap_dynamics(data = longitudinalMetabolomics)

## ----fig.aling='center',fig.dpi=150-------------------------------------------
longitudinalMetabolomics <- compare_metabolites(
  data = longitudinalMetabolomics
)
heatmap_metabolites(data = longitudinalMetabolomics)

## ----fig.aling='center',fig.dpi=150-------------------------------------------
dynamics <- metadata(longitudinalMetabolomics)[["comparison_dynamics"]][["estimates"]]
metabolites <- metadata(longitudinalMetabolomics)[["comparison_metabolites"]]
temp <- left_join(dynamics, metabolites, join_by("cluster_a", "cluster_b"))
x <- unique(c(temp[, "cluster_a"], temp[, "cluster_b"]))
temp <- temp %>% mutate(scale_Jaccard = scale(Jaccard))
ggplot(temp, aes(x = cluster_b, y = cluster_a)) +
  geom_point(aes(size = Jaccard, col = mu_mean)) +
  theme_bw() +
  scale_color_viridis_c(option = "magma") +
  scale_x_discrete(limits = x) +
  xlab("") +
  ylab("") +
  scale_y_discrete(limits = x) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  labs(col = "dynamics distance", size = "metabolite similarity") +
  ggtitle("comparison of clusters")

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