Skip to content

GitLab

Commit d5ae43b3 authored by Marek Ostaszewski's avatar Marek Ostaszewski
Browse files

Cleaned crosstalks code, Reactome added, new Cytoscape session file

parent d9709d11
Pipeline #37217 passed with stages
in 3 minutes and 17 seconds
Hide whitespace changes
Inline Side-by-side
Showing with 229 additions and 157 deletions
Resources/Crosstalks/find_crosstalks.R
View file @ d5ae43b3
......@@ -7,16 +7,46 @@
### Helper functions to access and retrieve content from MINERVA-hosted diagrams
source("https://git-r3lab.uni.lu/covid/models/-/raw/master/Integration/MINERVA_access/minerva_access_functions.R")
### For handling WikiPathway queries
library(rWikiPathways)
### A function retrieving element identifiers from MINERVA and WikiPathway diagrams
get_grouped_refs <- function(fcomponents, fgroups = NULL) {
### For handling Reactome pathways
### Get HGNCs of UniProts
translate_up_hgnc <- function(upid, hard_filter = T) {
resp <- ask_GET(paste0("https://www.ebi.ac.uk/proteins/api/genecentric/", upid))
if(is.null(resp)) { return("") }
res <- fromJSON(resp)
if(is.null(res$relatedGene)) { return("") }
### This is a patchy reference; the correct way would be to get HGNC ID (numeric) and retrieve the symbol based on this
### but I'd like to avoid yet another API call
hgncs <- with(res$relatedGene, geneName[accession == upid & geneNameType == "Gene name"])
### Gene symbols retrieved this way may not be HGNCs; we hard filter out non-uppercased or too short ids
if(hard_filter) {
hgncs <- hgncs[(nchar(hgncs) > 2) & (toupper(hgncs) == hgncs)]
}
return(hgncs)
}
### Get participants of a pathway and translate them into HGNCs
get_reactome_pathway_contents <- function(reactome_path_id) {
res <- fromJSON(ask_GET(paste0("https://reactome.org/ContentService/data/participants/", reactome_path_id)))
ids <- sapply(res$refEntities, "[[", "identifier")
symbols <- sapply(unique(unlist(ids)), translate_up_hgnc)
return(unlist(unique(symbols)))
}
### A function retrieving element identifiers from MINERVA
get_mnv_refs <- function(fcomponents) {
### Get MINERVA HGNCs as a named list
message("Retrieving MINERVA HGNCs collection")
refs <- sapply(fcomponents$map_elements, function(x) get_annotation(x$references, "HGNC_SYMBOL"))
names(refs) <- fcomponents$models$name
### Get WikiPathways HGNCs as a named list
mnv_refs <- sapply(fcomponents$map_elements, function(x) get_annotation(x$references, "HGNC_SYMBOL"))
names(mnv_refs) <- fcomponents$models$name
return(mnv_refs)
}
### A function retrieving element identifiers from WikiPathway diagrams
get_wp_refs <- function() {
message("Retrieving WP collection")
### This file contains all the pathways in the WP COVID-19 repo
wpids <- readLines("https://raw.githubusercontent.com/wikipathways/SARS-CoV-2-WikiPathways/master/pathways.txt")
......@@ -24,12 +54,39 @@ get_grouped_refs <- function(fcomponents, fgroups = NULL) {
message("Retrieving WP HGNCs")
wp_hgncs <- sapply(wpids, getXrefList, "H")
names(wp_hgncs) <- wpids
refs <- c(refs, wp_hgncs)
### If no grouping needed, return the lists as they are
if(is.null(fgroups)) { return(refs) }
### If by-name groupings provided, group unique HGNCs and return
grefs <- sapply(fgroups, function(x) setdiff(unique(unlist(refs[x])), c("", NA)))
names(grefs) <- names(fgroups)
return(wp_hgncs)
}
### A function retrieving element identifiers from Reactome diagrams
get_rct_refs <- function() {
message("Retrieving Reactome collection")
rct_cov_hgncs <- sapply(c("R-HSA-9678110", "R-HSA-9679504", "R-HSA-9679514", "R-HSA-9683701", "R-HSA-9679509"),
get_reactome_pathway_contents)
rct_cov2_hgncs <- sapply(c("R-HSA-9694614", "R-HSA-9694676", "R-HSA-9694682", "R-HSA-9694635", "R-HSA-9694322"),
get_reactome_pathway_contents)
return(c(rct_cov_hgncs, rct_cov2_hgncs))
}
### Beautification of the graph output: WikiPathway names
name_wp_nicely <- function(wp_ids) {
### For each WP id, grab the "name" parameter from the output of the getPathwayInfo()
sapply(wp_ids, function(x)
paste0(x, ": ", rWikiPathways::getPathwayInfo(x)$name
))
}
### Beautification of the graph output: Reactome names
name_rct_nicely <- function(rct_ids) {
### For each Reactome id, grab the "displayName" parameter from the JSON describing the pathway
sapply(rct_ids, function(x)
paste0(x, ": ",
fromJSON(ask_GET(paste0("https://reactome.org/ContentService/data/query/enhanced/",x)))$displayName))
}
### A function grouping the references according to a given structure
group_refs <- function(refs, groups) {
### Using by-name groupings provided, group unique HGNCs and return
grefs <- sapply(groups, function(x) setdiff(unique(unlist(refs[x])), c("", NA)))
names(grefs) <- names(groups)
return(grefs)
}
......@@ -92,14 +149,19 @@ get_external_crosstalks <- function(frefs, external_interactions) {
### A function operating on the external cross-talks table (see above)
### following the format: from, to, from pathway, from pathway n, to pathway, to pathway n
### list external edges to and from pathways
get_external_edges <- function(fcrosstalk) {
get_external_edges <- function(fcrosstalk, simplify = T) {
ext_edgelist <- c()
for(i in 1:nrow(fcrosstalk)) {
from <- unlist(strsplit(fcrosstalk[i,3],";"))
for(j in from) { ext_edgelist <- rbind(ext_edgelist, c(j, fcrosstalk[i,1])) }
ext_edgelist <- rbind(ext_edgelist, c(fcrosstalk[i,1], fcrosstalk[i,2]))
to <- unlist(strsplit(fcrosstalk[i,5],";"))
for(j in to) { ext_edgelist <- rbind(ext_edgelist, c(fcrosstalk[i,2], j)) }
if(simplify) { ### Same direction: pathway to element
for(j in from) { ext_edgelist <- rbind(ext_edgelist, c(j, fcrosstalk[i,1])) }
for(j in to) { ext_edgelist <- rbind(ext_edgelist, c(j, fcrosstalk[i,2])) }
} else {
for(j in from) { ext_edgelist <- rbind(ext_edgelist, c(j, fcrosstalk[i,1])) }
for(j in to) { ext_edgelist <- rbind(ext_edgelist, c(fcrosstalk[i,2], j)) }
}
}
return(unique(ext_edgelist))
}
......@@ -109,47 +171,21 @@ map <- "https://covid19map.elixir-luxembourg.org/minerva/api/"
### Get the components of MINERVA and WikiPathways (see minerva_access.R)
map_components <- get_map_components(map)
### Introduce pathway groups (following 'https://covid.pages.uni.lu/map_contents')
map_groups <- list(
`Virus replication cycle` = c(
"WP4846", "WP4799",
"Virus replication cycle", "Nsp9 protein interactions", "SARS-CoV-2 RTC and transcription",
"Nsp4 and Nsp6 protein interactions","E protein interactions", "Orf3a protein interactions"),
`ER stress and unfolded protein response` = c(
"WP4861",
"Endoplasmatic Reticulum stress", "Orf10 Cul2 pathway"),
`Autophagy and protein degradation` = c(
"WP4860", "WP4936", "WP4863",
"Endoplasmatic Reticulum stress"),
`Apoptosis pathway` = c(
"WP4864", "WP4877",
"Apoptosis pathway", "JNK pathway", "Electron Transport Chain disruption"),
`Renin-angiotensin system` = c(
"WP4883", "WP4965", "WP4969",
"Renin-angiotensin pathway"),
`Coagulopathy pathway` = c(
"WP4927",
"Coagulation pathway"),
`PAMP signalling` = c(
"WP4912",
"PAMP signalling"),
`Induction of interferons and the cytokine storm` = c(
"WP4868", "WP4880", "WP4876", "WP4884", "WP4961",
"Interferon 1 pathway", "Interferon lambda pathway"),
`Altered host metabolism` = c(
"WP4853", "WP4904",
"Kynurenine synthesis pathway", "HMOX1 pathway", "Nsp14 and metabolism", "Pyrimidine deprivation"),
`Other diagrams` = c(
"WP4891", "WP5017", "TGFbeta signalling", "overview"
)
)
map_groups <- fromJSON("map_groups.json")
### Detailed lists
mnv_refs <- get_mnv_refs(map_components)
wp_refs <- get_wp_refs()
rct_refs <- get_rct_refs()
### For a non-grouped version, re-run without 'map_groups'
map_grefs <- get_grouped_refs(map_components, map_groups)
### Get existing between-pathway crosstalks
ix_edgelist <- get_crosstalk_edges(map_grefs, compress = F)
### Global list of all identifiers
all_hgncs <- unique(unlist(map_grefs))
### Grouped list
map_grefs <- group_refs(c(mnv_refs, wp_refs, rct_refs), map_groups)
### For downstream display of graphs
names(mnv_refs) <- paste0("MINERVA: ", names(mnv_refs))
names(wp_refs) <- name_wp_nicely(names(wp_refs))
names(rct_refs) <- name_rct_nicely(names(rct_refs))
map_refs <- c(mnv_refs, wp_refs, rct_refs)
###
### For new interactions, we will use OmniPathDB (https://omnipathdb.org) and
......@@ -206,26 +242,39 @@ emmaa_interactors <- t(unname(sapply(emmaa_hgnc, function(x) c(x$subj$name, x$ob
### Add filtering by the OmniPathDB contents
emmaa_op <- paste0(emmaa_interactors[,1],"-",emmaa_interactors[,2]) %in%
paste0(op_narrow$source_genesymbol,"-", op_narrow$target_genesymbol)
emmaa_interactors_op <- emmaa_interactors[emmaa_op,]
### New crosstalks based on EMMAA TM results, raw and filtered by OmniPathDB
ix_emmaa <- get_external_crosstalks(map_grefs, emmaa_interactors)
ix_emmaa_op <- get_external_crosstalks(map_grefs, emmaa_interactors[emmaa_op,])
### Create a network of new inter-pathway crosstalks (with more than one 'from_pathway_n' and 'to_pathway_n')
### (raw and OmniPathDB-filtered)
ix_emmaa_edges <- get_external_edges(ix_emmaa[ix_emmaa[,'from_pathway_n'] > 1 &
ix_emmaa[,'to_pathway_n'] > 1,])
ix_emmaa_edges_op <- get_external_edges(ix_emmaa_op[ix_emmaa_op[,'from_pathway_n'] > 1 &
ix_emmaa_op[,'to_pathway_n'] > 1,])
### Create a network of new (HGNC not among the map contents) upstream regulators (with 'from_pathway_n' == 0, and 'to_pathway_n' > 0)
### (raw and OmniPathDB-filtered)
ix_reg_emmaa_edges <- get_external_edges(ix_emmaa[ix_emmaa[,'from_pathway_n'] == 0 &
ix_emmaa[,'to_pathway_n'] > 0 &
!(ix_emmaa[,'from'] %in% all_hgncs),])
ix_reg_emmaa_edges_op <- get_external_edges(ix_emmaa_op[ix_emmaa_op[,'from_pathway_n'] == 0 &
ix_emmaa_op[,'to_pathway_n'] > 0 &
!(ix_emmaa_op[,'from'] %in% all_hgncs),])
### Generate results for different graph combinations
### selected_refs are either pathways or groups
### selected_emmaa_interactors can be NULL (for existing crosstalks), raw (without op) or filtered (with op)
### new_elements T for new regulators or F for new crosstalks
generate_combination <- function(selected_refs, selected_emmaa_interactors = NULL, new_elements = F) {
if(is.null(selected_emmaa_interactors)) {
### Get existing crosstalks between selected references (pathways/groups)
return(get_crosstalk_edges(selected_refs, compress = F))
} else {
### Combine pathway/group interactions with selected EMMAA interactions
ix_refs_emmaa <- get_external_crosstalks(selected_refs, selected_emmaa_interactors)
if(new_elements) {
### New upstream regulators: return new interactions that have
### 1. no incoming interactions from existing elements
### 2. at least one outgoing interaction to an existing element
### 3. are not among existing elements
return(get_external_edges(ix_refs_emmaa[ix_refs_emmaa[,'from_pathway_n'] == 0 &
ix_refs_emmaa[,'to_pathway_n'] > 0 &
!(ix_refs_emmaa[,'from'] %in% unique(unlist(selected_refs))),]))
} else {
### Existing crosstalks: return new interactions that have
### 1. more than one incoming interaction from an existing element
### 2. more than one outgoing interaction to an existing element
return(get_external_edges(ix_refs_emmaa[ix_refs_emmaa[,'from_pathway_n'] > 1 &
ix_refs_emmaa[,'to_pathway_n'] > 1,]))
}
}
return(NULL)
}
### Use igraph to create xgmml files (for CytoScape) and pdf files
library(igraph)
......@@ -235,23 +284,16 @@ library(igraph)
produce_graph <- function(fedgelist, fxgmml = NULL, fpdf = NULL) {
### Create a graph from a non-empty names
f_g <- graph.edgelist(fedgelist[fedgelist[,1] != "",])
### Graphical properties
### Pathway nodes contain spaces in their names, we need them for a different visualisation
pathway_indices <- grep(" ", V(f_g)$name)
### Basic HGNC node colour
f_color <- rep("lightblue", gorder(f_g))
### Three-neighbour colour
if(quantile(degree(f_g))[2] > 2) { f_color[degree(f_g) >= quantile(degree(f_g))[2]] <- "yellow" }
### Four-or-more-neighbour colour
if(quantile(degree(f_g))[3] > 3) { f_color[degree(f_g) >= quantile(degree(f_g))[3]] <- "red" }
### Pathway colour
f_color[pathway_indices] <- "lightgreen"
### Height (for the rectangle shape)
f_height <- rep(40, gorder(f_g))
f_color <- rep("lightblue", gorder(f_g)) ### Basic HGNC node colour
f_color[degree(f_g) >= 3] <- "yellow" ### Three-neighbour colour
f_color[degree(f_g) >= 5] <- "red" ### Five-or-more-neighbour colour
f_color[grep(" ", V(f_g)$name)] <- "lightgreen" ### Pathway colour
f_height <- rep(35, gorder(f_g)) ### Height (for the rectangle shape)
f_height[grep(" ", V(f_g)$name)] <- 45 ### Pathway colour
f_g <- set.vertex.attribute(f_g, "color", value = f_color) ### Set the colours
### Set the colours
f_g <- set.vertex.attribute(f_g, "color", value = f_color)
### Create layout, stretch a bit
l = layout_with_fr(f_g)
l[,1] <- l[,1]*2000
......@@ -268,15 +310,26 @@ produce_graph <- function(fedgelist, fxgmml = NULL, fpdf = NULL) {
return(f_g)
}
### Create a graph of direct crosstalks
ix_g <- produce_graph(ix_edgelist, fxgmml = "ix_out.graphml", fpdf = "ix_out.pdf")
### Create a graph of new EMMAA crosstalks (raw)
ix_emmaa_g <- produce_graph(ix_emmaa_edges, fxgmml = "ix_emmaa_out.graphml", fpdf = "ix_emmaa_out.pdf")
### Create a graph of new EMMAA crosstalks (OmniPathDB-filtered)
ix_emmaa_op_g <- produce_graph(ix_emmaa_edges_op, fxgmml = "ix_emmaa_op_out.graphml", fpdf = "ix_emmaa_op_out.pdf")
### Create a graph of new EMMAA upstream regulators (raw)
ix_reg_emmaa_g <- produce_graph(ix_reg_emmaa_edges, fxgmml = "ix_reg_emmaa_out.graphml", fpdf = "ix_reg_emmaa_out.pdf")
### Create a graph of new EMMAA upstream regulators (OmniPathDB-filtered)
ix_reg_emmaa_op_g <- produce_graph(ix_reg_emmaa_edges_op, fxgmml = "ix_reg_emmaa_op_out.graphml", fpdf = "ix_reg_emmaa_op_out.pdf")
### Here, only a selection of graphs is built
### (other combinations can be created by using "generate_combination()" with other parameters)
### 1. Pathway existing crosstalks
ix_g <- produce_graph(generate_combination(selected_refs = map_refs),
fxgmml = "ix.graphml", fpdf = "ix.pdf")
### 2. Group existing crosstalks
grix_g <- produce_graph(generate_combination(selected_refs = map_grefs),
fxgmml = "grix.graphml", fpdf = "grix.pdf")
### 3. Pathway new crosstalks EMMAA + OmniGraffle
ix_emmaa_op_g <- produce_graph(generate_combination(selected_refs = map_refs,
selected_emmaa_interactors = emmaa_interactors_op),
fxgmml = "ix_emmaa_op.graphml", fpdf = "ix_emmaa_op.pdf")
### 4. Group new crosstalks EMMAA + OmniGraffle
grix_emmaa_op_g <- produce_graph(generate_combination(selected_refs = map_grefs,
selected_emmaa_interactors = emmaa_interactors_op),
fxgmml = "grix_emmaa_op.graphml", fpdf = "grix_emmaa_op.pdf")
### 5. Group new upstream regulators EMMAA + OmniGraffle
greg_emmaa_op_g <- produce_graph(generate_combination(selected_refs = map_grefs,
selected_emmaa_interactors = emmaa_interactors_op,
new_elements = T),
fxgmml = "greg_emmaa_op.graphml", fpdf = "greg_emmaa_op.pdf")
message("Done.")
\ No newline at end of file
Resources/Crosstalks/map_groups.json 0 → 100644
View file @ d5ae43b3
{
"Virus replication cycle": [
"WP4846",
"WP4799",
"Virus replication cycle",
"Nsp9 protein interactions",
"SARS-CoV-2 RTC and transcription",
"Nsp4 and Nsp6 protein interactions",
"E protein interactions",
"Orf3a protein interactions",
"R-HSA-9678110",
"R-HSA-9679504",
"R-HSA-9694514",
"R-HSA-9683701",
"R-HSA-9679509",
"R-HSA-9694614",
"R-HSA-9694676",
"R-HSA-9694682",
"R-HSA-9694635",
"R-HSA-9694322"
],
"ER stress and unfolded protein response": [
"WP4861",
"WP5027",
"Endoplasmatic Reticulum stress",
"Orf10 Cul2 pathway"
],
"Autophagy and protein degradation": [
"WP4860",
"WP4936",
"WP4863",
"Endoplasmatic Reticulum stress"
],
"Apoptosis pathway": [
"WP4864",
"WP4877",
"WP5038",
"Apoptosis pathway",
"JNK pathway",
"Electron Transport Chain disruption"
],
"Renin-angiotensin system": [
"WP4883",
"WP4965",
"WP4969",
"Renin-angiotensin pathway"
],
"Coagulopathy pathway": [
"WP4927",
"Coagulation pathway"
],
"PAMP signalling": [
"WP4912",
"PAMP signalling"
],
"Induction of interferons and the cytokine storm": [
"WP4868",
"WP4880",
"WP4876",
"WP4884",
"WP4961",
"WP5039",
"Interferon 1 pathway",
"Interferon lambda pathway",
"TGFbeta signalling",
"NLRP3 inflammasome activation"
],
"Altered host metabolism": [
"WP4853",
"WP4904",
"Kynurenine synthesis pathway",
"HMOX1 pathway",
"Nsp14 and metabolism",
"Pyrimidine deprivation"
],
"Other diagrams": [
"WP4891",
"WP5017",
"overview"
]
}
Resources/Crosstalks/minerva_access.R deleted 100644 → 0
View file @ d9709d11
##################################################
## Project: COVID-19 Disease Map
## Script purpose: Convenience functions for accessing the MINERVA Platform
## Date: 24.12.2020
## Author: Marek Ostaszewski
##################################################
library(httr)
library(jsonlite)
### A convenience function to handle API queries
ask_GET <- function(fask_url, verbose = F) {
if(verbose) {
message(URLencode(fask_url))
}
resp <- httr::GET(url = URLencode(fask_url),
httr::add_headers('Content-Type' = "application/x-www-form-urlencoded"))
if(httr::status_code(resp) == 200) {
return(httr::content(resp, as = "text"))
}
return(NULL)
}
### Get the components of a given map/project on the MINERVA Platform
get_map_components <- function(map_api, project_id = NULL) {
if(is.null(project_id)) {
### If project id not given, get configuration of the map, to obtain the latest (default) version
cfg <- fromJSON(ask_GET(paste0(map_api, "configuration/")))
project_id <- cfg$options[cfg$options$type == "DEFAULT_MAP","value"]
}
### The address of the latest (default) build
mnv_base <- paste0(map_api, "projects/",project_id,"/")
message(paste0("Asking for diagrams in: ", mnv_base, "models/"))
### Get diagrams
models <- fromJSON(ask_GET(paste0(mnv_base, "models/")))
### Get elements of the chosen diagram
model_elements <- lapply(models$idObject,
function(x)
fromJSON(ask_GET(paste0(mnv_base,"models/",x,"/",
"bioEntities/elements/?columns=id,name,type,references,elementId,complexId")),
flatten = F))
### Request for reactions that have at least one top 10 element as participant
model_reactions <- lapply(models$idObject,
function(x)
fromJSON(ask_GET(paste0(mnv_base,"models/",x,"/",
"bioEntities/reactions/?columns=modifiers,products,reactants")),
flatten = F))
### Pack all into a list and return
return(list(models = models, map_elements = model_elements, map_reactions = model_reactions))
}
### Get annotation of a given type, from element/reaction references
get_annotation <- function(freferences, ftype) {
sapply(freferences,
function(x) {
ifelse(any(x$type == ftype), x$resource[x$type == ftype], NA)
})
}
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment