---
title: "Using MetaboDynamics with data frames"
package: MetaboDynamics
author: "Katja Danielzik (katja.danielzik@uni-due.de)"

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

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

# Background 
The MetaboDynamics vignette explains the package workflow and interpretation of results with a SummarizedExperiment object as input.
In this vignette I will show the usage of MetaboDynamics with a data frame as input. 

# Setup: load required packages
```{r setup}
library(MetaboDynamics)
library(SummarizedExperiment) # storing and manipulating simulated metabolomics data
library(ggplot2) # visualization
library(dplyr) # data handling
library(tidyr) # data handling
```

# Load data
The MetaboDynamics package also includes a simulated example data set in a data frame format. 
To load the data set in data frame format, execute the following code:

```{r}
data("longitudinalMetabolomics_df", package = "MetaboDynamics")
```

# Model dynamics
To keep run time of this vignette short we will execute the workflow on a subset
of five metabolites across all conditions. 

```{r,fig.wide=TRUE}
data("longitudinalMetabolomics_df")
# take a sample (of five metabolites and subset data
samples <- c(
  "UMP", "PEP",
  "2-Aminomuconate","ATP"
)
data <- longitudinalMetabolomics_df %>% filter(metabolite %in% samples)

head(data)
```
The measured metabolites are stored in column "metabolite", time points are 
specified in column "time", and column "condition" specifies the experimental condition.

The required log-transformed and scaled (per metabolite and condition to a mean 
of 0 and sd of 1) metabolite abundances are stored in column "m_scaled". The
function fit_dynamics_model() allows to specify other column names for metabolite,
time points, condition and scaled measurement. 

First we will fit the dynamics model and extract the diagnostics:

```{r}
# fit dynamics model
fit <- #  when using a data frame as input fits have to stored in a separate object
  fit_dynamics_model(
    model = "scaled_log",
    data = data,
    scaled_measurement = "m_scaled", # in which column are scaled measurments stored?
    max_treedepth = 10,
    adapt_delta = 0.99, # default 0.95
    iter = 2000,
    warmup = 2000 / 4, # default is 1/4 of iterations
    chains = 1, # only set to 1 in vignette, recommended default is 4!
    cores = 1 # only set to 1 in vignette, can be same as chains if machine allows for parallelization
  )
```

This returns a model fit.

```{r}
# Extract diagnostics
diagnostics <- # when using a data frame as input diagnostics have to be stored in a separate object
  diagnostics_dynamics(
    data = data, # data frame that was used to fit the dynamics model,
    fit = fit, # list of fits from dynamics model, result of fit_dynamics_mode function
    iter = 2000, # how many iterations were used to fit the dynamics model
    chains = 1, # how many chains were used to fit the dynamics model
  )
```

That returns a list with elements [["model_diagnostics"]] which holds a summary
of the model diagnostics and the extracted posterior estimates. [["posterior]]
holds the posterior predictions needed for the posterior predictive check (PPC)

To visualize the diagnostics and the PPC the following code can be used:
```{r}
plot_diagnostics(
  data = data, # data frame used to fit the dynamics model
  diagnostics = diagnostics[["model_diagnostics"]] # summary of diagnostics
)
plot_PPC(
  data = data, posterior = diagnostics[["posterior"]],
  scaled_measurement = "m_scaled"
)
```
In this case the diagnostics are all indicating a sane model fit: The model fitting did not result in 
divergent transitions, the maximum tree depth was not exceeded, Rhat values were
below 1.01 and number of effective samples exceeded 100. 

# Extract estimates and visualize results
To extract the estimates one has to specify again the data, list of fits and number of iterations and chains used to fit the model. Additionally is is possible to specify the columns in which metabolite names, experimental conditions and time points are stored. This is equivalent to the fit_dynamics_model function.

Additionally one can specify how many samples should be drawn from the posterior
to for example test clustering performance. 

```{r}
estimates <- # estimates have to be stored in a separate object when using data frames
  estimates_dynamics(
    data = data,
    fit = fit)
```

This returns a list containing the estimates for every experimental condition.
The estimates can be visualized in two ways by the following code:

```{r}
# only visualize differences between time points
plot_estimates(
  # does not need data as input
  estimates = estimates,
  delta_t = TRUE, # choose to visualize differences between time points
  dynamics = FALSE,
  distance_conditions = FALSE
)

# only visualize dynamics
plot_estimates(
  estimates = estimates,
  delta_t = FALSE,
  distance_conditions = FALSE,
  dynamics = TRUE
) # choose to visualize the dynamics

# only visualize euclidean distance between dynamics vectors
# only visualize dynamics
plot_estimates(
  estimates = estimates,
  delta_t = FALSE,
  distance_conditions = TRUE,
  dynamics = FALSE
) # choose to visualize the dynamics
```

# Clustering dynamics

To cluster the dynamics we need either the estimates from the dynamics model (stored in object "estimates") 
or a data frame in which the means of measurements are stored in a column named "mu_mean", categorical time points are stored in ascending order in column "time.ID" and experimental conditions in a column named "condition".

The following chunk of code shows the needed data frame format for the clustering function:

```{r}
head(estimates)
```

With the output from the estimates_dynamics function we can cluster 
metabolite dynamics per experimental condition. 

```{r}
cluster <- # clustering results have to be stored in separate object when using data frame as input
  cluster_dynamics(
    estimates = estimates, # data is now the estimates or a data frame ob similar structure
    fit = fit,
    distance = "euclidean", # which distance method should be used
    agglomeration = "ward.D2", # which agglomeration method for hierarchical clustering should be used
    deepSplit = 2, # sensitivity of cluster analysis,
    minClusterSize = 1, # minimum number of metabolites in one cluster
    B = 10, # number of bootstrapps
    )
```

To visualize the clustering results we can use the function plot_cluster which
returns a list of plots. For details see Vignette "MetaboDynamics".

```{r}
plots <- plot_cluster(cluster)
plots$trees$A
plots$clusterplots$A
plots$lineplots$A
plots$patchwork$A
```


# Over-Representation Analysis (ORA)
To conduct the over-representation analysis KEGG IDs of the experimental metabolites are needed.
Additionally, information on all metabolites that are stored in the KEGG database is needed. This KEGG background information is stored in the data set "modules_compounds". To apply ORA to a specific dataset a second dataset "metabolite_modules" is needed that can be obtained by filtering "modules_compounds" for the experimental metabolites.

The clustering result also needs to be converted from a list to a single data frame.

```{r}
# load background dataset
data("modules_compounds")
data("metabolite_modules")
# get KEGGs from data
data("IDs")
```

The over-representation analysis and visualization can be achieved with the
following code: 

```{r}
ora <- # store ORA result to separate object when using data frames as input
  ORA_hypergeometric(
    background = modules_compounds,
    annotations = metabolite_modules,
    data = cluster, # dataframe format of clustering output
    tested_column = "middle_hierarchy",
    IDs = IDs
  )

# Visualization
plot_ORA(
  data = ora,
  tested_column = "middle_hierarchy"
)
```

# Compare dynamics clusters
For the comparison of dynamics clusters between experimental conditions one
can employ two apporaches: comparing dynamics between clusters and comparing
metabolite composition between clusters.

## Dynamics
To compare dynamics between clusters awhile using a data frame as input use the following
code:

```{r}
comparison_dynamics <- # result needs to be stored in a separate object when using data frames
  compare_dynamics(
    data = cluster,
    # clustering result
    cores = 1 # only set to 1 for vignette, can be increased to 4 for parallelization
  )
```

The data frame needed for visualization of the results is the list element [["estimates"]]
of the function results. To visualize the results run the following code:

```{r}
# Visualize comparison results
heatmap_dynamics(
  estimates = comparison_dynamics[["estimates"]],
  data = cluster
)
```


## Metabolites
The comparison of metabolite composition follows the same principle as the comparison of dynamics between clusters:

```{r}
# compare metabolite composition
compare_metabolites <-
  compare_metabolites(
    data = cluster
  )

# Visualization
heatmap_metabolites(
  distances = compare_metabolites,
  data = cluster
)
```


## Combine both
The combination of both comparisons may facilitate detection of differences of longitudinal metabolomes between experimental conditions. 

```{r}
# combine comparison results
temp <- left_join(
  comparison_dynamics[["estimates"]], # dynamics comparison
  compare_metabolites,
  join_by("cluster_a", "cluster_b") # join by cluster comparisons
)

# get unique clusters
x <- unique(c(temp[, "cluster_a"], temp[, "cluster_b"]))

# draw plot
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")
```


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