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.
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")
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:
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: #>  viper: matrix with 2505 rows and 50 columns #>  cnv: matrix with 598 rows and 50 columns #>  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.683252
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] 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) #>  "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-06
Initialize the Moma object with the assays and pathway data in order to start the analysis. The required inputs are the following:
Other optional inputs include:
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 matrix
The first step,
runDIGGIT() will run the DIGGIT inference algorithm [paper] to
find statistical interactions between VIPER-inferred proteins and genomic events.
makeInteractions() function will infer robust computational predictions
using all the provided data, including the Conditional Inference of Network Dynamics
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`$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
The results of this analysis can be accessed directly in the following result fields:
These will be used for plotting the genomic saturation curves as well.
cluster1.checkpoint <- momaObj$checkpoints[] print (cluster1.checkpoint[1:5]) #>  "5702" "55170" "7014" "23028" "29128"
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.
The other important results of the analysis are the statisically significant
genomic events found to be upstream of each of the checkpoint master regulators.
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:
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[])
plots$curve.plots[] + ggtitle("Genomic Saturation Curve for GBM Subtype 3")
Any of the results fields can be saved by using the
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")
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: #>  LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C #>  LC_TIME=en_US.UTF-8 LC_COLLATE=C #>  LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 #>  LC_PAPER=en_US.UTF-8 LC_NAME=C #>  LC_ADDRESS=C LC_TELEPHONE=C #>  LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C #> #> attached base packages: #>  grid stats graphics grDevices utils datasets methods #>  base #> #> other attached packages: #>  ggplot2_3.3.1 MOMA_1.0.2 BiocStyle_2.16.0 #> #> loaded via a namespace (and not attached): #>  Biobase_2.48.0 tidyr_1.1.0 #>  splines_4.0.0 BiocManager_1.30.10 #>  stats4_4.0.0 GenomeInfoDbData_1.2.3 #>  yaml_2.2.1 robustbase_0.93-6 #>  pillar_1.4.4 lattice_0.20-41 #>  glue_1.4.1 limma_3.44.1 #>  digest_0.6.25 GenomicRanges_1.40.0 #>  RColorBrewer_1.1-2 XVector_0.28.0 #>  qvalue_2.20.0 colorspace_1.4-1 #>  htmltools_0.4.0 Matrix_1.2-18 #>  plyr_1.8.6 pkgconfig_2.0.3 #>  GetoptLong_0.1.8 magick_2.3 #>  bookdown_0.19 zlibbioc_1.34.0 #>  purrr_0.3.4 scales_1.1.1 #>  MultiAssayExperiment_1.14.0 tibble_3.0.1 #>  farver_2.0.3 generics_0.0.2 #>  IRanges_2.22.2 ellipsis_0.3.1 #>  withr_2.2.0 SummarizedExperiment_1.18.1 #>  BiocGenerics_0.34.0 magrittr_1.5 #>  crayon_1.3.4 evaluate_0.14 #>  tools_4.0.0 MKmisc_1.6 #>  hms_0.5.3 GlobalOptions_0.1.1 #>  lifecycle_0.2.0 matrixStats_0.56.0 #>  ComplexHeatmap_2.4.2 stringr_1.4.0 #>  S4Vectors_0.26.1 munsell_0.5.0 #>  cluster_2.1.0 DelayedArray_0.14.0 #>  compiler_4.0.0 GenomeInfoDb_1.24.0 #>  rlang_0.4.6 RCurl_1.98-1.2 #>  rjson_0.2.20 circlize_0.4.9 #>  labeling_0.3 bitops_1.0-6 #>  rmarkdown_2.2 gtable_0.3.0 #>  reshape2_1.4.4 R6_2.4.1 #>  knitr_1.28 dplyr_1.0.0 #>  clue_0.3-57 shape_1.4.4 #>  readr_1.3.1 stringi_1.4.6 #>  parallel_4.0.0 Rcpp_126.96.36.199 #>  vctrs_0.3.0 png_0.1-7 #>  DEoptimR_1.0-8 tidyselect_1.1.0 #>  xfun_0.14