OmnipathR 3.2.8
1 Institute for Computational Biomedicine, Heidelberg University
In this tutorial we use Gene Symbols (HGNC symbols) as primary identifiers just to maintain better human readability across all intermediate steps, however the same procedure can be done with UniProt IDs, and this latter is more preferable.
The endpoints of pathways are transcription factors (TFs). We obtain a list of TFs from Vaquerizas et al. 2009, the so called TF census:
library(OmnipathR)
library(magrittr)##
## Attaching package: 'magrittr'## The following object is masked from 'package:tidyr':
##
## extractlibrary(dplyr)
tfcensus <-
import_omnipath_annotations(
resources = 'TFcensus',
entity_types = 'protein'
) %>%
pull(genesymbol) %>%
uniqueAlternatively we can directly download the same dataset, but retrieving data from OmniPath is preferable most of the times:
library(purrr)##
## Attaching package: 'purrr'## The following object is masked from 'package:magrittr':
##
## set_names## The following objects are masked from 'package:igraph':
##
## compose, simplifytfcensus_direct <-
tfcensus_download() %>%
pull(`HGNC symbol`) %>%
discard(is.na) %>%
uniqueAnother way to obtain a list of TFs is to collect the proteins with known regulons (target genes). We can download the TF-target interactions from OmniPath and collect all proteins which are source of any of these interactions.
tfs_from_network <-
import_transcriptional_interactions(
# confidence levels;
# we use only the 3 levels with highest confidence
dorothea_levels = c('A', 'B', 'C'),
entity_types = 'protein'
) %>%
pull(source_genesymbol) %>%
uniqueThe latest method resulted only around 600 TFs, however 97 of these are not in TF census. Maybe a good solution is to use the union of these and TF census:
tfs <- union(tfcensus, tfs_from_network)OmniPath contains a few resources providing subcellular localization data. We present here a three of them:
localizations <-
import_omnipath_annotations(
resources = c(
'Ramilowski_location',
'UniProt_location',
'HPA_subcellular'
),
entity_types = 'protein',
wide = TRUE
)Checking how many proteins are covered in each of these resources:
localizations %>%
map(function(x){x %>% pull(genesymbol) %>% n_distinct})## $HPA_subcellular
## [1] 6802
##
## $Ramilowski_location
## [1] 18866
##
## $UniProt_location
## [1] 16765The Intercell database of OmniPath is a comprehensive resource of cell-cell communication roles. We obtain a set of ligands, receptors, and as alternatives we create a sets of transmitters and receivers, which are broader categories, generalizations of the former two. After loading each of these datasets, we print the number of unique proteins.
ligands <-
OmnipathR::import_omnipath_intercell(
parent = 'ligand',
topology = 'sec',
consensus_percentile = 50,
loc_consensus_percentile = 30,
entity_types = 'protein'
) %>%
pull(genesymbol) %>%
unique %T>%
length
receptors <-
OmnipathR::import_omnipath_intercell(
parent = 'receptor',
topology = 'pmtm',
consensus_percentile = 50,
loc_consensus_percentile = 30,
entity_types = 'protein'
) %>%
pull(genesymbol) %>%
unique %T>%
length
transmitters <-
OmnipathR::import_omnipath_intercell(
causality = 'trans',
topology = 'sec',
consensus_percentile = 50,
loc_consensus_percentile = 30,
entity_types = 'protein'
) %>%
pull(genesymbol) %>%
unique %T>%
length
receivers <-
OmnipathR::import_omnipath_intercell(
causality = 'rec',
topology = 'pmtm',
consensus_percentile = 50,
loc_consensus_percentile = 30,
entity_types = 'protein'
) %>%
pull(genesymbol) %>%
unique %T>%
lengthThe protein-protein interactions in OmniPath come from many resources and
are divided into 4 major datasets (omnipath, pathwayextra, kinaseextra
and ligrecextra). We can try our procedure with any combination of these,
or even select only some resources within them. Here, for the sake of
simplicity, we use only the core omnipath PPI dataset, which contains
only literature curated, causal interactions.
ppi <-
import_omnipath_interactions(
datasets = 'omnipath',
entity_types = 'protein'
)Many databases contains curated canonical pathways, such as SIGNOR, KEGG or SignaLink. SignaLink pathways have been curated by experts from ligands to TFs, and the proteins are annotated as ligands, receptors, mediators, scaffolds, co-factors and TFs. For this reason, we use SignaLink pathways in this demonstration.
signalink_pathways <-
import_omnipath_annotations(
resources = 'SignaLink_pathway',
entity_types = 'protein',
wide = TRUE
)
signalink_functions <-
import_omnipath_annotations(
resources = 'SignaLink_function',
entity_types = 'protein',
wide = TRUE
)In this example, we use the ligands of the SignaLink TGF-beta pathway as genes of interest, which will be the starting point of the paths:
ligands_of_interest <-
intersect(
signalink_pathways %>%
filter(pathway == 'TGF') %>%
pull(genesymbol),
signalink_functions %>%
filter(`function` == 'Ligand') %>%
pull(genesymbol)
) %T>%
lengthThe data frame of interactions has to be converted to an igraph graph
object, because the find_all_paths function works on this kind of
object.
ppi_graph <-
ppi %>%
interaction_graphSearching paths in the graph, we have to define the source and target nodes and the maximum length of the paths. The source nodes will be the ligands of the TGF-beta pathway, the target nodes the TFs, and the maximum path length I set to 3. This latter results paths of 3 edges, connecting 4 nodes. The more numerous are the source and target node sets, the longer is the maximum length, the longer it takes to look up all the paths.
First we make sure we use only the ligands and TFs which are part of this network:
library(igraph)
ligands_of_interest__in_network <-
ligands_of_interest %>%
keep(. %in% V(ppi_graph)$name)
tfs__in_network <-
tfs %>%
keep(. %in% V(ppi_graph)$name)Then we look up the paths, this might take a few minutes:
paths <-
ppi_graph %>%
find_all_paths(
start = ligands_of_interest__in_network,
end = tfs__in_network,
maxlen = 3,
attr = 'name'
)A maxlen of 2 resulted around 600 paths, increasing the length only by 1
results 15,173 paths and takes ca. 10 min to run. A length of 4 probably
would be still feasible, but might result millions of paths and hours of
computation time.
All the resulted paths start at a ligand and end on a TF, but we need to make sure that the intermediate nodes meet the desired constraints: the second node is a plasma membrane receptor, the mediators are located in the cytoplasm, etc.
# only selecting Cytoplasm, of course we could
# consider other intracellular compartments
in_cytoplasm <-
localizations$UniProt_location %>%
filter(location == 'Cytoplasm') %>%
pull(genesymbol) %>%
unique
in_nucleus <-
localizations$UniProt_location %>%
filter(location %in% c('Nucleus', 'Nucleolus')) %>%
pull(genesymbol) %>%
unique
paths_selected <-
paths %>%
# removing single node paths
discard(
function(p){
length(p) == 1
}
) %>%
# receptors are plasma membrane transmembrane
# according to our query to OmniPath
keep(
function(p){
p[2] %in% receptors
}
) %>%
# making sure all further mediators are
# in the cytoplasm
keep(
function(p){
all(p %>% tail(-2) %>% is_in(in_cytoplasm))
}
) %>%
# the last nodes are all TFs, implying that they are in the nucleus
# but we can optionally filter foall(p %>% tail(-2)r this too:
keep(
function(p){
last(p) %in% in_nucleus
}
) %>%
# finally, we can remove paths which contain TFs also as intermediate
# nodes as these are redundant
discard(
function(p){
!any(p %>% head(-1) %>% is_in(tfs))
}
)Out of the 15,173 paths, 2,045 (~13%) meet these criteria. 43 of the paths consist of 3 nodes (ligand->receptor->TF), while 2,002 of 4 nodes (ligand->receptor->mediator->TF):
paths_selected %>% length
paths_selected %>% map_int(length) %>% tableWe have seen find_all_paths is limited in efficiency.
A more convenient method could be to use only data frames and dplyr’s
join methods to build the paths. This comes with the advantage that
even more efficient backends are available for dplyr, such as dtplyr
(data.frame backend) or dbplyr (relational DBMS backends), hence if
needed the performance can be easily improved.
In OmnipathR we have already a dedicated
function for creating ligand-receptor (or transmitter-receiver) interactions,
we can use this as a first step. For simplicity, here we limited the query to
ligand-receptor interactions and do a simple filtering to select high
confidence ones (mostly based on consensus across resources). Alternatively,
one can fine tune better the parameters to import_intercell_network, or
call it without any parameters, and use filter_intercell_network for
quality filtering. At the end we use simplify_intercell_network to get
rid of redundant columns.
ligand_receptor <-
import_intercell_network(
ligand_receptor = TRUE,
high_confidence = TRUE,
entity_types = 'protein'
) %>%
simplify_intercell_networkWe query the TF-target gene interactions into another data frame:
tf_target <-
import_transcriptional_interactions(
# confidence levels;
# we use only the 2 levels with highest confidence
dorothea_levels = c('A', 'B'),
# I added this only to have less interactions so the
# examples here run faster; feel free to remove it,
# then you will have more gene regulatory interactions
# from a variety of databases
datasets = 'tf_target',
entity_types = 'protein',
# A workaround to mitigate a temporary issue (05/01/2022)
resources = c('ORegAnno', 'PAZAR')
)The PPI network we already have in the ppi data frame, we loaded it in the
previous section. If you missed it, you can load it like this:
ppi <-
import_omnipath_interactions(
datasets = 'omnipath',
entity_types = 'protein'
)As a next step, we add UniProt localization annotations to the network data frames, except the ligand-receptor data frame, as this is created in a way that all ligands are extracellular and all receptors are plasma membrane transmembrane. In addition, we also add kinase annotations from the kinase.com resource, just to showcase working with multiple annotations. One protein might have multiple localizations and belong to multiple kinase families. We transform the localization and kinase family columns to list columns, which means one field might contain more than one values.
ppi %<>%
annotated_network('UniProt_location') %>%
annotated_network('kinase.com') %>%
# these columns define a unique interaction; 5 of them would be enough
# for only the grouping, we include the rest only to keep them after
# the summarize statement
group_by(
source,
target,
source_genesymbol,
target_genesymbol,
is_directed,
is_stimulation,
is_inhibition,
sources,
references,
curation_effort,
n_references,
n_resources
) %>%
summarize(
location_source = list(unique(location_source)),
location_target = list(unique(location_target)),
family_source = list(unique(family_source)),
family_target = list(unique(family_target)),
subfamily_source = list(unique(subfamily_source)),
subfamily_target = list(unique(subfamily_target))
) %>%
ungroup## `summarise()` has grouped output by 'source', 'target', 'source_genesymbol', 'target_genesymbol',
## 'is_directed', 'is_stimulation', 'is_inhibition', 'sources', 'references', 'curation_effort', 'n_references'.
## You can override using the `.groups` argument.tf_target %<>%
annotated_network('UniProt_location') %>%
annotated_network('kinase.com') %>%
group_by(
source,
target,
source_genesymbol,
target_genesymbol,
is_directed,
is_stimulation,
is_inhibition,
sources,
references,
curation_effort,
# Workaround to mitigate a temporary issue with
# DoRothEA data in OmniPath (05/01/2022)
# dorothea_level,
n_references,
n_resources
) %>%
summarize(
location_source = list(unique(location_source)),
location_target = list(unique(location_target)),
family_source = list(unique(family_source)),
family_target = list(unique(family_target)),
subfamily_source = list(unique(subfamily_source)),
subfamily_target = list(unique(subfamily_target))
) %>%
ungroup## `summarise()` has grouped output by 'source', 'target', 'source_genesymbol', 'target_genesymbol',
## 'is_directed', 'is_stimulation', 'is_inhibition', 'sources', 'references', 'curation_effort', 'n_references'.
## You can override using the `.groups` argument.We build the paths starting from the ligand-receptor interactions, joining
the PPI interactions repeatedly until we reach the desired path length.
After each step we throw away the paths which contain loops; set aside the
paths which end in a TF; and keep growing the rest of the paths. At the end,
we add the gene regulatory interactions to the end of each path, making all
paths start in an extracellular ligand and end in a TF target gene. In each
step we add suffixes to the column names: _step1, _step2, and so on, to
ensure the attributes of each step can be recognized. Each step will extend
the width of the data frame by around 20 columns. We could reduce this by
dropping certain variables. First, we define the set of TFs:
tfs <-
tf_target %>%
pull(source_genesymbol) %>%
uniqueFor the repetitive task of growing paths we define a function:
library(rlang)##
## Attaching package: 'rlang'## The following objects are masked from 'package:purrr':
##
## %@%, as_function, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl, flatten_raw,
## invoke, splice## The following object is masked from 'package:magrittr':
##
## set_names## The following object is masked from 'package:igraph':
##
## is_namedlibrary(tidyselect)
grow_paths <- function(paths, step, interactions, endpoints = NULL){
right_target_col <- sprintf(
'target_genesymbol_step%i',
step
)
paths$growing %<>% one_step(step, interactions)
# finished paths with their last node being an endpoint:
paths[[sprintf('step%i', step)]] <-
paths$growing %>%
{`if`(
is.null(endpoints),
.,
filter(., !!sym(right_target_col) %in% endpoints)
)}
# paths to be further extended:
paths$growing %<>%
{`if`(
is.null(endpoints),
NULL,
filter(., !(!!sym(right_target_col) %in% endpoints))
)}
invisible(paths)
}
one_step <- function(paths, step, interactions){
left_col <- sprintf('target_genesymbol_step%i', step - 1)
right_col <- sprintf('source_genesymbol_step%i', step)
by <- setNames(right_col, left_col)
paths %>%
# making sure even at the first step we have the suffix `_step1`
{`if`(
'target_genesymbol' %in% colnames(.),
rename_with(
.,
function(x){sprintf('%s_step%i', x, step)}
),
.
)} %>%
{`if`(
step == 0,
.,
inner_join(
.,
# adding suffix for the next step
interactions %>%
rename_with(function(x){sprintf('%s_step%i', x, step)}),
by = by
)
)} %>%
# removing loops
filter(
select(., contains('_genesymbol')) %>%
pmap_lgl(function(...){!any(duplicated(c(...)))})
)
}This function seems to be useful, after some refactoring I will add it to
OmnipathR. But now let’s use it for building the paths. First, using
reduce, we call the previously defined function repeatedly according to
the desired number of steps. After, we go through the collections of paths
of various lengths, and add the last step, the TF-target interactions.
This procedure is extremely memory intensive, if you run it with 3 steps,
make sure you have at least 12G of free memory. Higher number of steps
(longer paths) require even more memory. Of course this can be mitigated
by removing some columns, adding less annotations, or using smaller networks
(more stringent queries or subsequent filtering).
library(stringr)
steps <- 2 # to avoid OOM
paths <-
seq(0, steps) %>%
head(-1) %>%
reduce(
grow_paths,
interactions = ppi,
endpoints = tfs,
.init = list(
growing = ligand_receptor
)
) %>%
within(rm(growing)) %>%
map2(
names(.) %>% str_extract('\\d+$') %>% as.integer %>% add(1),
one_step,
interactions = tf_target
)Here is a list of all packages used in this vignette:
library(OmnipathR)
library(magrittr)
library(dplyr)
library(purrr)
library(igraph)
library(rlang)
library(tidyselect)
library(stringr)
library(rmarkdown)This way you can make sure all these packages are installed:
missing_packages <-
setdiff(
c(
'magrittr',
'dplyr',
'purrr',
'igraph',
'rlang',
'tidyselect',
'stringr',
'rmarkdown',
'devtools'
),
installed.packages()
)
if(length(missing_packages)){
install.packages(missing_packages)
}
if(!'OmnipathR' %in% installed.packages()){
library(devtools)
devtools::install_github('saezlab/OmnipathR')
}This is an Rmd (R markdown) document, a popular reporting format for R.
You can generate various document formats from it, such as HTML or PDF,
using the render function from the rmarkdown package:
library(rmarkdown)
render('paths.Rmd', output_format = 'html_document')Below you find information about the operating system, R version and all the loaded packages and their versions from an R session where this vignette run without errors. This can be useful for troubleshooting, as different software versions are one possible cause of errors.
sessionInfo()## R version 4.1.2 (2021-11-01)
## 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 LC_TIME=en_GB
## [4] LC_COLLATE=C LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] stringr_1.4.0 tidyselect_1.1.2 rlang_1.0.1 purrr_0.3.4 magrittr_2.0.2 gprofiler2_0.2.1
## [7] dnet_1.1.7 supraHex_1.32.0 hexbin_1.28.2 tidyr_1.2.0 ggraph_2.0.5 igraph_1.2.11
## [13] ggplot2_3.3.5 dplyr_1.0.8 OmnipathR_3.2.8 BiocStyle_2.22.0
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-155 bitops_1.0-7 bit64_4.0.5 progress_1.2.2 httr_1.4.2
## [6] Rgraphviz_2.38.0 tools_4.1.2 backports_1.4.1 bslib_0.3.1 utf8_1.2.2
## [11] R6_2.5.1 lazyeval_0.2.2 DBI_1.1.2 BiocGenerics_0.40.0 colorspace_2.0-3
## [16] withr_2.4.3 gridExtra_2.3 prettyunits_1.1.1 bit_4.0.4 curl_4.3.2
## [21] compiler_4.1.2 graph_1.72.0 cli_3.2.0 xml2_1.3.3 plotly_4.10.0
## [26] labeling_0.4.2 bookdown_0.24 sass_0.4.0 scales_1.1.1 checkmate_2.0.0
## [31] readr_2.1.2 rappdirs_0.3.3 digest_0.6.29 rmarkdown_2.11 pkgconfig_2.0.3
## [36] htmltools_0.5.2 fastmap_1.1.0 highr_0.9 htmlwidgets_1.5.4 readxl_1.3.1
## [41] jquerylib_0.1.4 farver_2.1.0 generics_0.1.2 jsonlite_1.8.0 vroom_1.5.7
## [46] RCurl_1.98-1.6 Matrix_1.4-0 Rcpp_1.0.8 munsell_0.5.0 fansi_1.0.2
## [51] ape_5.6-1 logger_0.2.2 viridis_0.6.2 lifecycle_1.0.1 stringi_1.7.6
## [56] yaml_2.3.5 MASS_7.3-55 grid_4.1.2 parallel_4.1.2 ggrepel_0.9.1
## [61] crayon_1.5.0 lattice_0.20-45 graphlayouts_0.8.0 hms_1.1.1 magick_2.7.3
## [66] knitr_1.37 pillar_1.7.0 stats4_4.1.2 glue_1.6.1 evaluate_0.15
## [71] data.table_1.14.2 BiocManager_1.30.16 vctrs_0.3.8 tzdb_0.2.0 tweenr_1.0.2
## [76] cellranger_1.1.0 gtable_0.3.0 polyclip_1.10-0 assertthat_0.2.1 xfun_0.29
## [81] ggforce_0.3.3 tidygraph_1.2.0 later_1.3.0 viridisLite_0.4.0 tibble_3.1.6
## [86] ellipsis_0.3.2