---
title: "MetaboDynamics: analyzing longitudinal metabolomics data with probabilistic models"
package: MetaboDynamics
author: "Katja Danielzik (katja.danielzik@uni-due.de)"

output:
  BiocStyle::html_document:
    toc: true
    toc_depth: 2
vignette: >
  %\VignetteIndexEntry{1. MetaboDynamics}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r style, echo = FALSE, results = 'asis'}
BiocStyle::markdown()
```

# Background 
This package was developed to facilitate the analysis of longitudinal 
metabolomics data. Most tools only allow the comparison between two time points 
or experimental conditions and are using frequentist statistical methods.

As metabolomics data is often noisy and we generally have few replicates due
to high costs, a robust method is needed for the estimation of mean
concentrations at every time point. For this we employ a Bayesian hierarchical 
model that assumes normal distributions of log-transformed metabolite 
concentrations.

Here we want to show a complete workflow to analyze concentration tables with 
probabilistic models.

# Installation
MetaboDynamics can be installed from the devel branch of Bioconductor. Open R
4.4 and execute the following code:

```{r,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: load required packages
```{r setup}
library(MetaboDynamics)
library(SummarizedExperiment)
library(ggplot2)
library(dplyr)
library(tidyr)
```

# Load data and plot data overview

We have a simulated data (called "longitudinalMetabolomics") set of 98 metabolites with 
three measurement replicates at four time points (1-4) across 3 experimental conditions (A-B).
In the first step in this workflow we estimate the dynamics of every single 
metabolite at every experimental condition (p.e. radiation dose).

The simulated data is represented as [SummarizedExperiment](https://bioconductor.org/packages/release/bioc/html/SummarizedExperiment.html) object where 
concentration tables of each experimental condition are stored in assays 
(raw concentrations in "concentrations", log-transformed transformations in 
"log_con" and scaled log-transformed concentrations" in "scaled_log") and 
metabolite names, KEGG IDs, experimental conditions and clustering solutions per
experimental condition are stored in colData and timepoint specifications in rowData.

The next plot shows the raw data.

```{r,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")
```

The raw data is not distributed normally. So let's log-transform the values.
In the integrated simulated data set this is already done in the column "log_m".

```{r,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")
```

The next plot shows the raw dynamics of single metabolites. 

```{r,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")
```

We define dynamics as deviations at the observed time points from the 
metabolite's mean concentration. As the raw concentrations of metabolites can 
differ by orders of magnitude from each other, and we want to be able to compare
dynamics of metabolites with each other, we standardize each metabolite at each
experimental condition to a mean of zero and a standard deviation of one. 

In the simulated data set the scaled measurements
are in column "m_scaled".

```{r,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")
```

Now we can finally model the dynamics. This might take of the order of 10 minutes 
per experimental condition for 98 metabolites. 

We employ a Bayesian hierarchical model with con = metabolite concentrations,
m = metabolite, c = experimental condition and t = time point ID:

\begin{align*}
\log(con_{m,c,t})&\sim {\sf normal}(\mu_{m,c,t},\sigma_{m,c,t}) \\ 
\mu_{m,c,t}&\sim {\sf normal}(0,2) \\
\sigma_{m,c,t}&\sim {\sf exponential}(\lambda_{m,c}) \\
\lambda_{m,c}&\sim {\sf exponential}(2) 
\end{align*}

The hierarchical Bayesian model provides a balance between over- and under fitting
of the model to the experimental data points.

The code below shows how to fit the model and how to extract the diagnostic
criteria from the model fits. 

# Model dynamics
```{r,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
)
```

This returns a list of model fits that are named by the experimental condition
("A","B","C"). If data is a summarizedExperiment the fits are stored in 
metadata(data)[["dynamic_fits"]]. 

With diagnostics_dynamics() we can extract all
the diagnostic criteria of MCMC runs to fit a Bayesian model (rhat, neff, 
divergences, max_treedepth) and visualize them. If data is a SummarizedExperiment
the diagnostics are stored in metadata(data)[["diagnostics_dynamics"]].
Additionally data frames for visual Posterior predictive checks (PPC) are 
prepared and plots for the PPCs and diagnostic criteria can be generated with 
plot_PPC() and plot_diagnostics().

```{r}
# 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
)
```

There should be no divergent transitions (number of divergent transitions=0)
and the maximum treedepth should not be exceeded (number of exceeded treedepth=0).

```{r}
# PPCs can be accessed with
plot_PPC(
  data = data
)
```

In the posterior predictive check all data points (black points) should lie
within the posterior of the Bayesian model (violin).

After checking the diagnostic criteria and the PPC we can extract the estimates:
If data is a SummarizedExperiment
the estimates are stored in metadata(data)[["estimates_dynamics"]]


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

We get two major outputs:
1. the estimation of concentration differences between two subsequent time points
of each metabolite at each experimental condition
2. the dynamic profiles of each metabolites at each experimental condition

## Differences between two timepoints
```{r,fig.wide=TRUE}
# 1) the differences between two timepoints
plot_estimates(
  data = data,
  dynamics = FALSE
)
```

If the 95% highest density interval (Credible Interval(CrI)) of the posterior 
does not include zero we can rather credibly state that there is a difference in 
mean concentrations between two time points. If the 95% HDI lies below zero we 
likely have a decrease in concentrations between the two time points, if it is 
above zero we likely have an increase in concentrations between time points. 

## Dynamic profiles
```{r,fig.wide=TRUE}
# 2) dynamic profiles
plot_estimates(
  data = data,
  delta_t = FALSE
)
```
So we now have dynamic profiles of many metabolites at each radiation dose.
We could now cluster these metabolite specific dynamics vectors 
(estimates[,c("mu1.mean":"mut.mean)]) to determine if groups of metabolites 
have similar dynamics.

# Cluster dynamics
For the sake of demonstration we use from here on a clustering result (data("cluster"))
on the full simulated data set (data("longitudinalMetabolomics")). In a real life example the 
optimal number of clusters ("k") should be determined by optimal clustering
criteria such as Gap statistics and average silhouette. 
The code below shows an example how the estimated dynamics profiles can be 
used for clustering. 

```{r,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
```

The metadata of the SummarizedExperiment object "longitudinalMetabolomics"
holds the clustering results of the whole simulated dataset data
("longitudinalMetabolomics").

```{r,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() +
  ylim(-2,2)+
  theme(legend.position = "none") +
  facet_grid(rows = vars(condition), cols = vars(cluster)) +
  ggtitle("clustered dynamics", "panels=cluster ID")
```

As we can see metabolites show different dynamics in different experimental
conditions. Can we quantify the biological function of these dynamics clusters?

# Over-representation analysis of functional modules in dynamics clusters
To quantify the possible biological function of these dynamics clusters we 
retrieved from the KEGG-database the following information with package KEGGREST:
1. to which functional modules our experimental metabolites are annotated and
2. which metabolites are annotated to functional modules in general.

The functional modules of the KEGG-database are organised in three hierarchies:
upper, middle and lower. Here we will conduct functional analysis on the middle 
hierarchy. To facilitate analysis the data frames "metabolite_modules", which 
holds the information about experimental metabolites, and "modules_compounds",
which holds the information about which metabolites are in general annotated to
functional modules, were prepared. The needed data frame for different sets
of experimental metabolites can be retrieved with the function get_ORA_annotations.
This function also allows to update the background KEGG information stored in 
"modules_compounds".

We load both data sets and can inspect the
documentation.

```{r}
data("metabolite_modules")
help("metabolite_modules")
head(metabolite_modules)
data("modules_compounds")
help("modules_compounds")
head(modules_compounds)
```

Here we have to keep in mind that not all KEGG modules are suitable for testing
on every observed organism and experimental condition. For example the modules
"Xenobiotics biodegradation","Biosynthesis of other secondary metabolites" and
"Biosynthesis of terpenoids and polyketides" should not be analyzed in a human
lung cancer cell line.

For the functional analysis we employ a hypergeometric model. We consider a 
functional module as over-represented in a cluster if the 95% inter-quantile
range (ICR) of the log-transformed probabilities of OvEs (observed vs expected) 
lies above zero. OvE refers to the ratio of observed metabolites in a cluster 
being mapped to a functional module over the number of expected metabolites in a
cluster being mapped to a module under the assumption of a hypergeometric 
distribution (=drawing without replacement). If data is a SummarizedExperiment
object where the clustering solution is stored in metadata(data)[["cluster"]]
the ORA results are stored in metadata(data)[["ORA_tested_column"]]
We apply the functional analysis to the middle and lower hierarchy of functional
modules. 

```{r,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")
```

Great, we can now see which functional module is over- 
(green points and error-bars) or under-represented (none in this example) in 
which dynamics cluster! For instance in cluster 3 at condition A and C the 
modules "Energy metabolism" and "Carbohydrate metabolism" are over-represented. 

# Comparison of clusters of different experimental conditions

## Dynamics 
We can not only conduct over-representation analysis of KEGG-functional modules but
also compare dynamics clusters across different experimental conditions. For 
this we employ a Bayesian model that estimates the mean difference as well as 
the standard deviation of differences between dynamics clusters. If data is a 
SummarizedExperiment object the results of compare_metabolites are stored 
in metadata(data)[["comparison_dynamics"]]

dist = vector of pairwise euclidean distances between each dynamic vector of 
cluster a and every dynamic vector of cluster b, ID = cluster pair ID
\begin{align*}
dist_{ID}&\sim {\sf normal}(\mu_{ID},\sigma_{ID}) \\
\mu_{ID}&\sim {\sf normal^+}(0,2) \\
\sigma_{ID}&\sim {\sf exponential}(1) 
\end{align*}

```{r,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
)
```

We can visualize the posterior estimates of the model with:

```{r,fig.aling='center',fig.dpi=150}
heatmap_dynamics(data = longitudinalMetabolomics)
```

The bigger and brighter a point, the smaller is the mean distance between
dynamics clusters and the smaller is the standard deviation. That means big 
bright points indicate high dynamic similarity which small spread. Here B_8 and
A_4 have high similarity in dynamics. 

## Metabolites
We can also compare metabolite composition of clusters. For this we employ
the Jaccard Index which is a metric for similarity of two vectors. Values near
0 indicate low similarity, values near 1 high similarity. If data is a 
SummarizedExperiment object the results of compare_metabolites are stored 
in metadata(data)[["comparison_metabolites"]].
```{r,fig.aling='center',fig.dpi=150}
longitudinalMetabolomics <- compare_metabolites(
  data = longitudinalMetabolomics
)
heatmap_metabolites(data = longitudinalMetabolomics)
```
The brighter a tile, the higher is the similarity of two clusters in regard
to their metabolite composition. 
We have two clusters that are similar in their metabolite composition:
C_6 and A_5. If we compare that to the dynamics profiles and ORA analysis
we see that similar functional modules are over-expressed as expected BUT
the dynamics differ between the two radiation doses. 

Can we facilitate visualization?

## Combine both
```{r,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")
```

We can find a cluster pair that is pretty similar in regards to their
composing metabolites but dissimilar in regards to their dynamics. Their ORA
profiles are quite similar as expected from the similar metabolite compositions 
but they show different dynamics between experimental conditions: B_7 and A_4

```{r}
sessionInfo()
```

