MOMA - Multi Omic Master Regulator Analysis
5 June 2020
Contents
1 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 // package.
2 Running MOMA
2.1 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.
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("MOMA")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)
example.gbm.mae
#> A MultiAssayExperiment object of 3 listed
#> experiments with user-defined names and respective classes.
#> Containing an ExperimentList class object of length 3:
#> [1] viper: matrix with 2505 rows and 50 columns
#> [2] cnv: matrix with 598 rows and 50 columns
#> [3] mut: matrix with 60 rows and 50 columns
#> Features:
#> experiments() - obtain the ExperimentList instance
#> colData() - the primary/phenotype DFrame
#> sampleMap() - the sample availability DFrame
#> `$`, `[`, `[[` - extract colData columns, subset, or experiment
#> *Format() - convert into a long or wide DFrame
#> assays() - convert ExperimentList to a SimpleList of matrices
MultiAssayExperiment::assays(example.gbm.mae)$viper[1:3, 1:3]
#> TCGA-14-1825-01 TCGA-16-0846-01 TCGA-32-2634-01
#> 165 -3.0369845 0.5837581 -1.982579
#> 196 -2.4667300 -3.7642532 -1.941096
#> 257 0.7872116 0.3260199 1.683252Next, 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] as well as protein-protein interactions predicted by the PrePPI structure-based algorithm [paper]. 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] 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.
names(gbm.pathways)
#> [1] "cindy" "preppi"
gbm.pathways$cindy[1:3]
#> $`55553`
#> 10000
#> 1
#>
#> $`171023`
#> 10013 58484
#> 5.257461e-06 1.000000e+00
#>
#> $`51295`
#> 10025 23370 51079 9546
#> 1.000000e+00 3.995671e-06 3.995671e-06 7.360446e-062.2 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]
- 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.
momaObj <- MomaConstructor(example.gbm.mae, gbm.pathways)
#> Found the following assays:viper, cnv, mut
#> Common samples across main assays: 50
#> Checking labels on pathway cindy
#> Found labels for 364 TFs in VIPER matrix
#> Checking labels on pathway preppi
#> Found labels for 2505 TFs in VIPER matrix2.3 MOMA Analysis on GBM Data
The first step, runDIGGIT() will run the DIGGIT inference algorithm [paper] 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.
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.
momaObj$Cluster()
#> using pearson correlation with weighted vipermat
#> testing clustering options, k = 2..15
#> Using reliability scores to select best solution.
#> Best solution is: 6clusters
# Explore the reliability scores
momaObj$clustering.results$all.cluster.reliability
#> 2clusters 3clusters 4clusters 5clusters 6clusters 7clusters 8clusters
#> 0.6144845 0.7529229 0.7456256 0.7775375 0.7907407 0.6766867 0.6616316
#> 9clusters 10clusters 11clusters 12clusters 13clusters 14clusters 15clusters
#> 0.6022523 0.5793794 0.5645946 0.5290190 0.5175175 0.5476476 0.5363463
# Save in the 3 cluster solution
momaObj$sample.clustering <- momaObj$clustering.results$`3clusters`$clusteringGenomic 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$checkpointsmomaObj$genomic.saturationmomaObj$coverage.summaryStats
These will be used for plotting the genomic saturation curves as well.
momaObj$saturationCalculation()cluster1.checkpoint <- momaObj$checkpoints[[1]]
print (cluster1.checkpoint[1:5])
#> [1] "5702" "55170" "7014" "23028" "29128"3 Plotting the Results
3.1 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.
3.2 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).
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.
library(ggplot2)
library(grid)
grid.draw(plots$oncoprint.plots[[3]])
plots$curve.plots[[3]] +
ggtitle("Genomic Saturation Curve for GBM Subtype 3")
4 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.
momaObj$saveData(output.folder = "outputs", "hypotheses", "checkpoints")5 Session Info
Here is the output of sessionInfo() on the system on which this package and documentation was compiled.
sessionInfo()
#> R version 4.0.0 (2020-04-24)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 18.04.4 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.11-bioc/R/lib/libRblas.so
#> LAPACK: /home/biocbuild/bbs-3.11-bioc/R/lib/libRlapack.so
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> attached base packages:
#> [1] grid stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] ggplot2_3.3.1 MOMA_1.0.2 BiocStyle_2.16.0
#>
#> loaded via a namespace (and not attached):
#> [1] Biobase_2.48.0 tidyr_1.1.0
#> [3] splines_4.0.0 BiocManager_1.30.10
#> [5] stats4_4.0.0 GenomeInfoDbData_1.2.3
#> [7] yaml_2.2.1 robustbase_0.93-6
#> [9] pillar_1.4.4 lattice_0.20-41
#> [11] glue_1.4.1 limma_3.44.1
#> [13] digest_0.6.25 GenomicRanges_1.40.0
#> [15] RColorBrewer_1.1-2 XVector_0.28.0
#> [17] qvalue_2.20.0 colorspace_1.4-1
#> [19] htmltools_0.4.0 Matrix_1.2-18
#> [21] plyr_1.8.6 pkgconfig_2.0.3
#> [23] GetoptLong_0.1.8 magick_2.3
#> [25] bookdown_0.19 zlibbioc_1.34.0
#> [27] purrr_0.3.4 scales_1.1.1
#> [29] MultiAssayExperiment_1.14.0 tibble_3.0.1
#> [31] farver_2.0.3 generics_0.0.2
#> [33] IRanges_2.22.2 ellipsis_0.3.1
#> [35] withr_2.2.0 SummarizedExperiment_1.18.1
#> [37] BiocGenerics_0.34.0 magrittr_1.5
#> [39] crayon_1.3.4 evaluate_0.14
#> [41] tools_4.0.0 MKmisc_1.6
#> [43] hms_0.5.3 GlobalOptions_0.1.1
#> [45] lifecycle_0.2.0 matrixStats_0.56.0
#> [47] ComplexHeatmap_2.4.2 stringr_1.4.0
#> [49] S4Vectors_0.26.1 munsell_0.5.0
#> [51] cluster_2.1.0 DelayedArray_0.14.0
#> [53] compiler_4.0.0 GenomeInfoDb_1.24.0
#> [55] rlang_0.4.6 RCurl_1.98-1.2
#> [57] rjson_0.2.20 circlize_0.4.9
#> [59] labeling_0.3 bitops_1.0-6
#> [61] rmarkdown_2.2 gtable_0.3.0
#> [63] reshape2_1.4.4 R6_2.4.1
#> [65] knitr_1.28 dplyr_1.0.0
#> [67] clue_0.3-57 shape_1.4.4
#> [69] readr_1.3.1 stringi_1.4.6
#> [71] parallel_4.0.0 Rcpp_1.0.4.6
#> [73] vctrs_0.3.0 png_0.1-7
#> [75] DEoptimR_1.0-8 tidyselect_1.1.0
#> [77] xfun_0.14