COSMOS-tutorial
A. Dugourd, A. Gabor and K. Zirngibl
11/10/2020
Installation
# 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
First, we load the package
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.
CARNIVAL_options <- cosmosR::default_CARNIVAL_options()
# CARNIVAL_options$solverPath <- "C:/Program Files/CPLEX_solver/cplex/bin/x64_win64/cplex.exe"
# CARNIVAL_options$solver <- "cplex" #or cbc
# 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 <- 2In the next section, we prepare the input to run cosmosR. The signaling inputs are the result of footprint based TF and kinase activity estimation. For more info on TF activity estimation 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_PUBCHEMidcompartment or XMetab_BIGGidcompartment or example “XMetab_6804m”. 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 it 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)
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
)## [1] "COSMOS: all 3 signaling nodes from data were found in the meta PKN"
## [1] "COSMOS: all 3 metabolic nodes from data were found in the meta PKN"
## [1] "COSMOS: 4660 of the 15919 genes in expression data were found as transcription factor target"
## [1] "COSMOS: 4660 of the 5312 transcription factor targets were found in expression data"
## [1] "COSMOS: removing unexpressed nodes from PKN..."
## [1] "COSMOS: 0 interactions removed"
## [1] "COSMOS: removing nodes that are not reachable from inputs within 15 steps"
## [1] "COSMOS: 86 from 115 interactions are removed from the PKN"
## [1] "COSMOS: 1 input/measured nodes are not in PKN any more: XMetab__439155___c____ and 0 more."
## [1] "COSMOS: removing nodes that are not observable by measurements within 15 steps"
## [1] "COSMOS: 10 from 29 interactions are removed from the PKN"
## [1] "COSMOS: 1 input/measured nodes are not in PKN any more: X1445 and 0 more."
## [1] "COSMOS: 0 interactions are removed from the PKN based on consistency check between TF activity and gene expression"
## [1] "lpSolve does not scale well with large PKNs. This solver is mainly for testing purposes. To run COSMSO, we recommend using cplex, or cbc solvers."
## [1] "lpSolve does not scale well with large PKNs. This solver is mainly for testing purposes. To run COSMSO, we recommend using cplex, or cbc solvers."
## [1] "COSMOS: 0 interactions are removed from the PKN based on consistency check between TF activity and gene expression"
## [1] "COSMOS: all 2 signaling nodes from data were found in the meta PKN"
## [1] "COSMOS: all 2 metabolic nodes from data were found in the meta PKN"
## [1] "COSMOS: 4660 of the 15919 genes in expression data were found as transcription factor target"
## [1] "COSMOS: 4660 of the 5312 transcription factor targets were found in expression data"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.
CARNIVAL_options$timelimit <- 14400
CARNIVAL_options$mipGAP <- 0.05
CARNIVAL_options$threads <- 2This is where cosmosR run.
test_result_for <- run_COSMOS_signaling_to_metabolism(data = test_for,
CARNIVAL_options = CARNIVAL_options)## [1] "lpSolve does not scale well with large PKNs. This solver is mainly for testing purposes. To run COSMSO, we recommend using cplex, or cbc solvers."
## [1] "lpSolve does not scale well with large PKNs. This solver is mainly for testing purposes. To run COSMSO, we recommend using cplex, or cbc solvers."Finally, we process the results of the first cosmosR run, to translate gene names and metabolites name.
data(metabolite_to_pubchem)
data(omnipath_ptm)
formated_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.
CARNIVAL_options$timelimit <- 3600
CARNIVAL_options$mipGAP <- 0.05
CARNIVAL_options$threads <- 2Now that the correct time is set up for the pre-optimisation run, we can prepare the inputs.
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)## [1] "COSMOS: all 3 signaling nodes from data were found in the meta PKN"
## [1] "COSMOS: all 3 metabolic nodes from data were found in the meta PKN"
## [1] "COSMOS: 4660 of the 15919 genes in expression data were found as transcription factor target"
## [1] "COSMOS: 4660 of the 5312 transcription factor targets were found in expression data"
## [1] "COSMOS: removing nodes that are not reachable from inputs within 15 steps"
## [1] "COSMOS: 105 from 115 interactions are removed from the PKN"
## [1] "COSMOS: 2 input/measured nodes are not in PKN any more: X4790, X5062 and 0 more."
## [1] "COSMOS: 2 input/measured nodes are not in PKN any more: XMetab__65359___c____, XMetab__107738___m____ and 0 more."
## [1] "COSMOS: removing nodes that are not observable by measurements within 15 steps"
## [1] "COSMOS: 5 from 10 interactions are removed from the PKN"
## [1] "COSMOS: 0 interactions are removed from the PKN based on consistency check between TF activity and gene expression"
## [1] "lpSolve does not scale well with large PKNs. This solver is mainly for testing purposes. To run COSMSO, we recommend using cplex, or cbc solvers."
## [1] "lpSolve does not scale well with large PKNs. This solver is mainly for testing purposes. To run COSMSO, we recommend using cplex, or cbc solvers."
## [1] "COSMOS: 0 interactions are removed from the PKN based on consistency check between TF activity and gene expression"
## [1] "COSMOS: all 1 signaling nodes from data were found in the meta PKN"
## [1] "COSMOS: all 1 metabolic nodes from data were found in the meta PKN"
## [1] "COSMOS: 4660 of the 15919 genes in expression data were found as transcription factor target"
## [1] "COSMOS: 4660 of the 5312 transcription factor targets were found in expression data"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.
## [1] "lpSolve does not scale well with large PKNs. This solver is mainly for testing purposes. To run COSMSO, we recommend using cplex, or cbc solvers."
## [1] "lpSolve does not scale well with large PKNs. This solver is mainly for testing purposes. To run COSMSO, we recommend using cplex, or cbc solvers."Finally we can format the result of the backward run as well (same as for forward run)
formated_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 overlap 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 overlap between the two solutions.
full_sif <- as.data.frame(rbind(formated_result_for[[1]], formated_result_back[[1]]))
full_attributes <- as.data.frame(rbind(formated_result_for[[2]], formated_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).
network_plot <- display_node_neighboorhood(central_node = 'NFKB1',
sif = full_sif,
att = full_attributes,
n = 7)
network_plotThis network represents 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), which in turn can be converted to glutathione GGT enzymes.
It is important to understand that each of this links is 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.
Tutorial section: Over Representation Analysis
Often it is useful to perform an Over Representation Analysis (ORA) on the resulting nodes of a COSMOS network as a first analysis step to get a more functional interpretation on the modeled signaling cascade. A common way to this is to test whether the selected genes (nodes) in the COSMOS solution network show statistically significant differences in comparison to the prior-knowledge network (PKN).
The differentially expressed genes give information about the cellular processes that are deregulated and if the proportions in various pathways are SIGNIFICANTLY different from what is expected.In this way the significant differences between two biological conditions (e.g. cancer vs. normal tissue, or treatment vs. untreated cells) can be shown.
Algorithms that perform an ORA are implemented in other R packages like piano or decoupleR. In addition to a gene set collection these algorithms require two different lists as inputs: - nodes in the COSMOS solution network which relate back to the input data (e.g. transcriptomics, proteomics, metabolomics, fluxomics, or perturbations) - all nodes (kinases, transcription factors, metabolites) in the prior-knowledge network (which are used as the background in our analysis)
In this section we will show how to obtain these two lists from a formated COSMOS result object.
sif_forward = formated_result_for[[1]]
att_forward = formated_result_for[[2]]
nodes_ORA = extract_nodes_for_ORA(
sif = sif_forward,
att = att_forward)Now this forground and background set can be used with any ORA anaylsis.
sessionInfo()## R version 4.1.1 (2021-08-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.14-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.14-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB 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] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] cosmosR_1.2.0
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.7 png_0.1-7 visNetwork_2.1.0
## [4] Biostrings_2.62.0 assertthat_0.2.1 digest_0.6.28
## [7] utf8_1.2.2 R6_2.5.1 GenomeInfoDb_1.30.0
## [10] stats4_4.1.1 RSQLite_2.2.8 evaluate_0.14
## [13] httr_1.4.2 pillar_1.6.4 zlibbioc_1.40.0
## [16] rlang_0.4.12 jquerylib_0.1.4 blob_1.2.2
## [19] S4Vectors_0.32.0 rmarkdown_2.11 readr_2.0.2
## [22] stringr_1.4.0 bcellViper_1.29.0 htmlwidgets_1.5.4
## [25] igraph_1.2.7 RCurl_1.98-1.5 bit_4.0.4
## [28] compiler_4.1.1 xfun_0.27 pkgconfig_2.0.3
## [31] BiocGenerics_0.40.0 CARNIVAL_2.4.0 htmltools_0.5.2
## [34] tidyselect_1.1.1 KEGGREST_1.34.0 tibble_3.1.5
## [37] GenomeInfoDbData_1.2.7 lpSolve_5.6.15 IRanges_2.28.0
## [40] fansi_0.5.0 crayon_1.4.1 dplyr_1.0.7
## [43] tzdb_0.2.0 bitops_1.0-7 jsonlite_1.7.2
## [46] lifecycle_1.0.1 DBI_1.1.1 magrittr_2.0.1
## [49] cli_3.1.0 stringi_1.7.5 vroom_1.5.5
## [52] cachem_1.0.6 XVector_0.34.0 bslib_0.3.1
## [55] ellipsis_0.3.2 generics_0.1.1 vctrs_0.3.8
## [58] dorothea_1.5.3 org.Hs.eg.db_3.14.0 tools_4.1.1
## [61] bit64_4.0.5 Biobase_2.54.0 glue_1.4.2
## [64] purrr_0.3.4 hms_1.1.1 parallel_4.1.1
## [67] fastmap_1.1.0 yaml_2.2.1 AnnotationDbi_1.56.0
## [70] memoise_2.0.0 knitr_1.36 sass_0.4.0