---
title: "MOMA - Multi Omic Master Regulator Analysis"
author: "Sunny Jones, Evan Paull, Mariano Alvarez, Federico Giorgi and Andrea Califano"
output: BiocStyle::html_document
#bibliography: bibliography.bib
vignette: >
  %\VignetteIndexEntry{MOMA - Multi Omic Master Regulator Analysis}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

# Introduction and Background

MOMA is a tool for inferring connections between Master Regulator proteins and 
genomic driver events in cancer. Master regulators are regulatory proteins 
(primarily transcription factors and co-transcription factors) that control cell
state. In the case of cancer and other disease states transcription factors have 
been shown to be key drivers of maintaining the disease state and can be targets 
for interventions. Often these master regulators are not mutated themselves but 
are downstream of mutations and other genomic alterations that ultimately 
dysregulate the normal activity of that regulator. MOMA uses multiple inputs of 
information to infer these connections and to improve the predictive value of 
the master regulator analysis.  



*For more information on Master Regulators and our tool for calculating their values see the paper published on VIPER in Nature Genetics in 2016* 
[paper](https://www.nature.com/articles/ng.3593) // [package](https://bioconductor.org/packages/release/bioc/html/viper.html).

*** 

# Running MOMA

## Getting Started

First install and load the library into the R session. Make sure all the other 
dependent packages have already been installed to ensure full functionality, 
including graphics and plotting functions.  


```{r setup, message = FALSE, results="hide", warning=FALSE, eval=FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("MOMA")
```


```{r load}
library(MOMA)
```


Explore the test data saved in the `MultiAssayExperiment` object. Confirm that we 
have 50 test samples and 2505 VIPER inferred proteins in the test viper matrix. 
This was generated by using `VIPER` on the GBM gene expression signature. The 
matrix has the samples across the columns and the regulators in the rows.  

The other assays needed for the experiment are:

- **SNP matrix**: binary 0/1 values where 1 indicates the prescence of a mutation
- **CNV matrix**: copy number variant values
- **Fusion matrix**: binary 0/1 values (optional)



```{r explore data, message = FALSE}
example.gbm.mae

MultiAssayExperiment::assays(example.gbm.mae)$viper[1:3, 1:3]
```

Next, we select the biological pathways we'd like to use to aid in our inference 
task. In this case, we're using the CINDy modulator inference algorithm [[paper](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0109569)] 
as well as protein-protein interactions predicted by the PrePPI structure-based 
algorithm [[paper](https://academic.oup.com/nar/article/41/D1/D828/1070064)]. 
The output of the CINDy algorithm is a likelihood (p value) for
the association of an upstream modulator with a particular regulators activity.
The relevant output for PrePPI is a likelihood (p value) that a modulator 
structurally binds to a regulator.

To make the MOMA object they need to be a named list of lists, with the indexes
of each being a regulator and their associated partners.

*note: CINDy values are context specific and need to be obtained for a particular tissue/expression set. Values calculated from the tumor types in TCGA can be found
[[here](https://figshare.com/articles/CINDy_Networks_from_TCGA/12198510)] for use in applying this analysis to other datasets. The PREPPi values are not context specific
and the values from this example analysis can be used with other datasets.* 

```{r pathways}
names(gbm.pathways)

gbm.pathways$cindy[1:3]
```

## Generating the **MOMA** object

Initialize the Moma object with the assays and pathway data in order to start
the analysis. The required inputs are the following:

- **Assays**: a MultiAssayExperiment object with the previously mentioned assays.
Alternatively a named list with the assays as matrices can be supplied. In both 
cased the names of the assays must be "viper", "cnv", "mut" and "fusion"
- **Pathway list**: see above

Other optional inputs include:

- **Gene Blacklist**: genes to not include in the analysis. By default the Broad
Institute's list of highly mutated genes are excluded [[link](https://software.broadinstitute.org/cancer/cga/mutsig)]
- **Output Folder** : Directory to save results to.
- **Gene Map** : Dataframe mapping Gene Names to Entrez IDs and Cytoband locations.
One is already provided by default. 

```{r moma object}
momaObj <- MomaConstructor(example.gbm.mae, gbm.pathways)
```

## MOMA Analysis on GBM Data

The first step, `runDIGGIT()`  will run the DIGGIT inference algorithm [[paper](https://www.cell.com/cell/fulltext/S0092-8674(14)01170-2)] to
find statistical interactions between VIPER-inferred proteins and genomic events.  

The `makeInteractions()` function will infer robust computational predictions
using all the provided data, including the Conditional Inference of Network Dynamics 
(CINDy) algorithm.  

The `Rank()` function will create a final ranking of candidate Master Regulators
for this cohort of patient samples.

```{r interactions, results = "hide", message = FALSE}
momaObj$runDIGGIT()

momaObj$makeInteractions()

momaObj$Rank()
```

Clustering of the samples, using the protein ranks computed in the last step, 
can then be performed using `Cluster()`. Multiple cluster solutions will be 
calculated, ranging from 2 to 15 clusters by default. The reliability or 
average silhouette scores of each can be assessed to determine an optimal 'k' 
number of clusters. By default the clustering solution with the maximum 
reliability will be saved to the object, but any solution can be saved in afterwards.

```{r clustering}
momaObj$Cluster()

# Explore the reliability scores
momaObj$clustering.results$all.cluster.reliability

# Save in the 3 cluster solution
momaObj$sample.clustering <- momaObj$clustering.results$`3clusters`$clustering

```

  
Genomic saturation analysis is then performed on each cluster with the 
`saturationCalculation()` function, allowing us to find the key proteins that are 
downstream of the majority of genomic events in the samples within a particular 
cluster. These regulators make up that cluster's `checkpoint`.  

The results of this analysis can be accessed directly in the following result 
fields:  
  
  - `momaObj$checkpoints` 
  - `momaObj$genomic.saturation`
  - `momaObj$coverage.summaryStats` 
  
These will be used for plotting the genomic saturation curves as well.
  
```{r saturation, message=FALSE, results="hide"}
momaObj$saturationCalculation()
```

```{r checkpoints}
cluster1.checkpoint <- momaObj$checkpoints[[1]]
print (cluster1.checkpoint[1:5])
```

***

# Plotting the Results

## Visualizing the VIPER matrix  

The primary results of the analysis are the master regululators of each particular
cluster's checkpoint, as displayed above. You can plot the original VIPER matrix 
subset down to only these regulators using whatever heatmap function you like. 


## Visualizing the Genomic Saturation and Events Plots

The other important results of the analysis are the statisically significant
genomic events found to be upstream of each of the checkpoint master regulators. 
The `makeSaturationPlots()` function takes the subtype specific genomic event 
interactions calculated in the previous step and makes plots for viewing the 
results. Because there are usually quite a large number of events that are 
detected only the most frequently occurring events will be plotted. 
(Note: to account for amplifications and deletions being events that can occur
across multiple genes, these events are considered on a cytoband basis).

```{r make plots, warning=FALSE, include=FALSE}
plots <- makeSaturationPlots(momaObj)
```


The two plot types are: 

- **Genomic Event Oncoprint Plots** : these plot the different types of 
genomic events found to be upstream of the regulators in the checkpoint,
using the style of cBioPortal's plotting. The samples in that cluster are in
the columns and the rows are the top significant events to be plotted (max
number can be changed). Summary bar plots are on the top and right side.
- **Saturation Curves** : these plot the number of regulators needed to capture
a certain percentage of genomic events. The inflection points plotted on
this graph match up with the cut off for which regulators make up a particular
cluster's checkpoint.   

The data for these plots are stored as grobs and ggplot objects so layers or other 
graphic modifications can be added afterwards for further customization.

```{r visualize plots}
library(ggplot2)
library(grid)
grid.draw(plots$oncoprint.plots[[3]])
plots$curve.plots[[3]] +
  ggtitle("Genomic Saturation Curve for GBM Subtype 3")
  
```

# Saving Data

Any of the results fields can be saved by using the `saveData()` function. 
Pass in any of the names of the results fields and they will be saved to files
in the designated output folder. If no names are passed all results will be saved.

```{r saving data, eval=FALSE}
momaObj$saveData(output.folder = "outputs", "hypotheses", "checkpoints")
```


# Session Info 

Here is the output of sessionInfo() on the system on which this package and
documentation was compiled.

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