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

## ----setup--------------------------------------------------------------------
library(MetaboDynamics)
library(SummarizedExperiment) # storing and manipulating simulated metabolomics data
library(ggplot2) # visualization
library(dplyr) # data handling
library(tidyr) # data handling

## -----------------------------------------------------------------------------
data("longitudinalMetabolomics", package = "MetaboDynamics")

## ----fig.wide=TRUE------------------------------------------------------------
# convert to dataframe
abundances <- as.data.frame(cbind(
  as.data.frame(t(assays(longitudinalMetabolomics)$concentrations)),
  as.data.frame(colData(longitudinalMetabolomics))
))
abundances <- abundances %>% pivot_longer(
  cols = seq_len(4), names_to = "time",
  values_to = "abundance"
)
ggplot(abundances, aes(x = abundance)) +
  geom_freqpoly() +
  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(abundances, aes(x = abundance)) +
  geom_density() +
  scale_x_log10() +
  theme_bw() +
  facet_grid(cols = vars(time), rows = vars(condition)) +
  ggtitle("data", "log-transformed values")

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

## ----fig.wide=TRUE------------------------------------------------------------
data("longitudinalMetabolomics")
# convert to dataframe
abundances <- as.data.frame(cbind(
  as.data.frame(t(assays(longitudinalMetabolomics)$scaled_log)),
  as.data.frame(colData(longitudinalMetabolomics))
))
abundances <- abundances %>% pivot_longer(
  cols = seq_len(4), names_to = "time",
  values_to = "abundance"
)
ggplot(abundances) +
  geom_line(aes(
    x = time,
    y = abundance, col = metabolite,
    group = interaction(metabolite, replicate)
  )) +
  theme_bw() +
  scale_color_viridis_d() +
  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
)

## -----------------------------------------------------------------------------
# get estimates
estimates <- metadata(data)[["estimates_dynamics"]]
# only show first rows of condition A
estimates[["A"]][1:10, ]

## ----fig.wide=TRUE------------------------------------------------------------
data <- cluster_dynamics(data,
                         deepSplit = 0) # low clustering sensitivity

## -----------------------------------------------------------------------------
plot <- plot_cluster(data)

## -----------------------------------------------------------------------------
plot[["A"]][["PCA_plot"]]

## -----------------------------------------------------------------------------
plot <- plot_cluster(longitudinalMetabolomics)
plot[["lineplot"]]

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

## ----fig.aling='center',fig.dpi=150-------------------------------------------
data("longitudinalMetabolomics")
longitudinalMetabolomics <- compare_dynamics(
  data = longitudinalMetabolomics,
  cores = 1 # just set to 1 for vignette, set to 4 if possible
)

## ----fig.aling='center',fig.dpi=150-------------------------------------------
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", "label = condition + cluster ID")

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