---
title: "COSMOS-tutorial"
author: "A. Dugourd, A. Gabor and K. Zirngibl"
date: "11/10/2020"
output:
  html_document: 
    toc: true
    toc_depth: 2
vignette: >
  %\VignetteIndexEntry{cosmosR tutorial}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
---
  
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Installation
```{r eval=FALSE}
# install from bioconductor
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("cosmosR")

# install the development version from GitHub
# install.packages("remotes")
remotes::install_github("saezlab/cosmosR")
```

## Introduction

COSMOS (Causal Oriented Search of Multi-Omic Space) is a method that integrates phosphoproteomics, transcriptomics, and metabolomics data sets. COSMOS leverages extensive prior knowledge of signaling pathways, metabolic networks, and gene regulation with computational methods to estimate activities of transcription factors and kinases as well as network-level causal reasoning. This pipeline can provide mechanistic explanations for experimental observations across multiple omic data sets.

![data_intro_figure](../inst/figures/intro_data.png)

First, we load the package

```{r, warning=FALSE, message=FALSE}
library(cosmosR)
```

## Tutorial section: signaling to metabolism

In this part, we can set up the options for the CARNIVAL run, such as timelimit and min gap tolerance.

The user should provide a path to its CPLEX/cbc executable

You can check the CARNIVAL_options variable to see all possible options that can be adjusted

In this example, we will use the built-in solver lpSolve. User should be aware that lpSolve should ONLY be used for TESTS. To obtain meaningful results, best solver is cplex, or cbc if not possible.

```{r, warning=FALSE, message=FALSE}
CARNIVAL_options <- cosmosR::default_CARNIVAL_options()
# CARNIVAL_options$solverPath <- "~/Documents/cplex"
# CARNIVAL_options$solver <- "cplex" #or cbc
CARNIVAL_options$solver <- "lpSolve" #or cbc
CARNIVAL_options$timelimit <- 3600
CARNIVAL_options$mipGAP <- 0.05
CARNIVAL_options$threads <- 2
```

In the next section, we prepare the input to run cosmosR
The signaling inputs are the result of footprint based TF and kinase activity estiamtion
For more info on TF activity estiamtion from transcriptomic data, see:https://github.com/saezlab/transcriptutorial (Especially chapter 4)

Here we use of toy PKN, to see the full meta PKN, you can load it with data(meta_network).

The metabolites in the prior knowledge network are identified as XMetab__PUBCHEMid___compartment____ or XMetab__BIGGid___compartment____
or example “XMetab__6804___m____”. The compartment code is the BIGG model standard (r, c, e, x, m, l, n, g). 
Thus we will first need to map whatever identifer for metabolite the data has to the one of the network.
Genes are identified as XENTREZid (in the signaling part of network) or XGene####__ENTREZid (in the reaction network part of network)

The maximum network depth will define the maximum number of step downstream of kinase/TF COSMOS will look for deregulated metabolites. Good first guess for max depth could be around 6 (here is 15 for the toy dataset)

The differential experession data is used to filter out wrong TF-target interactions in this context after a pre-optimisation.

The list of genes in the differential expression data will also be used as a reference to define which genes are expressed or not (all genes in the diff_expression_data are considered expressed, and genes that are no in diff_expression_data are removed from the network)

```{r, warning=FALSE, message=FALSE}
data(toy_network)
data(toy_signaling_input)
data(toy_metabolic_input)
data(toy_RNA)
test_for <- preprocess_COSMOS_signaling_to_metabolism(meta_network = toy_network,
                                        signaling_data = toy_signaling_input,
                                        metabolic_data = toy_metabolic_input,
                                                      diff_expression_data = toy_RNA,
                                                      maximum_network_depth = 15,
                                                      remove_unexpressed_nodes = TRUE,
                                                      CARNIVAL_options = CARNIVAL_options
                                                      )
```

In this part, we can set up the options for the actual run, such as timelimit and min gap tolerance.

The running time should be much higher here than in pre-optimisation. You cna increase the number of threads to use if you have many available CPUs.

```{r, warning=FALSE, message=FALSE}
CARNIVAL_options$timelimit <- 14400
CARNIVAL_options$mipGAP <- 0.05
CARNIVAL_options$threads <- 2
```

This is where cosmosR run.

```{r, warning=FALSE, message=FALSE}
test_result_for <- run_COSMOS_signaling_to_metabolism(data = test_for,
                                                      CARNIVAL_options = CARNIVAL_options)
```

Finally, we process the results of the first cosmosR run, to translate gene names and metabolites name.

```{r, warning=FALSE, message=FALSE}
data(metabolite_to_pubchem)
data(omnipath_ptm)
test_result_for <- format_COSMOS_res(test_result_for,
                                     metab_mapping = metabolite_to_pubchem,
                     measured_nodes = unique(c(names(toy_metabolic_input),
                                               names(toy_signaling_input))),
                                     omnipath_ptm = omnipath_ptm)
```

## Tutorial section: metabolism to signaling 

Before we run the metabolism to signaling part, we need to prepare again the inputs. 
```{r, warning=FALSE, message=FALSE}
CARNIVAL_options$timelimit <- 3600
CARNIVAL_options$mipGAP <- 0.05
CARNIVAL_options$threads <- 2
```

Now that the correct time is set up for the pre-optimisation run, we can prepare the inputs.

```{r, warning=FALSE, message=FALSE}
test_back <- preprocess_COSMOS_metabolism_to_signaling(meta_network = toy_network,
                                        signaling_data = toy_signaling_input,
                                        metabolic_data = toy_metabolic_input,
                                                       diff_expression_data = toy_RNA,
                                                       maximum_network_depth = 15,
                                                       remove_unexpressed_nodes = FALSE,
                                                       CARNIVAL_options = CARNIVAL_options
                                                       
)

```

Then we can run cosmosR to connect metabolism to signaling. The running time here usually needs to be longer, as this problem seems to be harder to solve for CPLEX.

```{r, warning=FALSE, echo=FALSE, message=FALSE}
CARNIVAL_options$timelimit <- 28800

test_result_back <- run_COSMOS_metabolism_to_signaling(data = test_back,
                                                       CARNIVAL_options = CARNIVAL_options)
```

Finally we can format the result of the backward run as well (same as for forward run)

```{r, warning=FALSE, message=FALSE}
test_result_back <- format_COSMOS_res(test_result_back,
                                      metab_mapping = metabolite_to_pubchem,
                      measured_nodes = unique(c(names(toy_metabolic_input),
                                                names(toy_signaling_input))),
                                      omnipath_ptm = omnipath_ptm)

```

## Tutorial section: Merge forward and backward networks and visualise network

Here we simply take the union of forward and backward runs to create a full network solution lopping between signaling, gene-regulation and metabolism. Since there is an overlapp between the result network of forward and backward run, you may optionally want  to check if there are any node sign that are incoherent in the overlapp between the two solutions.

```{r, warning=FALSE, message=FALSE}
full_sif <- as.data.frame(rbind(test_result_for[[1]], test_result_back[[1]]))
full_attributes <- as.data.frame(rbind(test_result_for[[2]], test_result_back[[2]]))

full_sif <- unique(full_sif)
full_attributes <- unique(full_attributes)
```

This function will generate a dynamic network plot centered on a given node of the network solution, and connecting it to measured nodes in the given range (here 5 steps).

```{r, warning=FALSE, message=FALSE}
network_plot <- display_node_neighboorhood(central_node = 'BCAT1', 
                                           sif = full_sif, 
                                           att = full_attributes, 
                                           n = 5)

network_plot
```
This network represent the flow of activities that can connect FOXM1 up-regulation with glutathione (CID 124886)
accumulation. Here, FOXM1 can activate MYC, which in turn activate BCAT1. The activation of BCAT1 can lead 
to the increased production of glutamate (CID 33032), whioch in turn can be converted to glutathione GGT enzymes.

It is important to understand that each of this links are hypothetical. The come from a larger pool of 
potential molecular interactions present in multiple online databases and compiled in omnipath, STITCH 
and recon metabolic network. They exist in the literature and are interactions that are known to
potentially exists in other experimental contexts. Thus, COSMOS compile all those potential interactions 
together and proposes a coherent set that can explain the data at hand. Here, such a set of mechanistic hypothesis
is plotted as a network connecting FOXM1 and glutathione production.

Those links should however be considered only as potential mechanistic connections, and will need to be further 
confirmed experimentally.
```{r, warning=FALSE, message=FALSE}
sessionInfo()
```
