Identifying the involvement in the same biological pathways
For each protein we retrieve the pathways they are involved in.
get_protein_pathway <- function(gene_protein_df) {
endpoint <- "/protein/in-pathway"
params <- list(
uniprot_id_list = gene_protein_df %>% pull(`uniprot_id`) %>% I()
)
df <- query_epigraphdb(route = endpoint, params = params, mode = "table", method = "POST")
if (nrow(df) > 0) {
res_df <- gene_protein_df %>%
select(`uniprot_id`) %>%
left_join(df, by = c("uniprot_id"))
} else {
res_df <- gene_protein_df %>%
select(`uniprot_id`) %>%
mutate(pathway_count = NA_integer_, pathway_reactome_id = NA_character_)
}
res_df <- res_df %>%
mutate(
pathway_count = ifelse(is.na(pathway_count), 0L, as.integer(pathway_count)),
pathway_reactome_id = ifelse(is.na(pathway_reactome_id), c(), pathway_reactome_id)
)
res_df
}
pathway_df <- get_protein_pathway(gene_protein_df = gene_protein_df)
pathway_df
#> # A tibble: 9 x 3
#> uniprot_id pathway_count pathway_reactome_id
#> <chr> <int> <list>
#> 1 Q13445 0 <NULL>
#> 2 Q9BYD3 0 <NULL>
#> 3 P29597 3 <chr [3]>
#> 4 Q8N9T8 0 <NULL>
#> 5 Q6P4F2 1 <chr [1]>
#> 6 O75791 3 <chr [3]>
#> 7 P05362 5 <chr [5]>
#> 8 Q9UMF0 2 <chr [2]>
#> 9 P0C6A0 0 <NULL>Now for each pair of proteins we match the pathways they have in common.
get_shared_pathway <- function(pathway_df) {
# For the protein-pathway data
# Get protein-protein permutations where they share pathways
per_permutation <- function(pathway_df, permutation) {
df <- pathway_df %>% filter(uniprot_id %in% permutation)
primary_pathway <- pathway_df %>%
filter(uniprot_id == permutation[1]) %>%
pull(pathway_reactome_id) %>%
unlist()
assoc_pathway <- pathway_df %>%
filter(uniprot_id == permutation[2]) %>%
pull(pathway_reactome_id) %>%
unlist()
intersect(primary_pathway, assoc_pathway)
}
pairwise_permutations <- pathway_df %>%
pull(`uniprot_id`) %>%
gtools::permutations(n = length(.), r = 2, v = .)
shared_pathway_df <- tibble(
protein = pairwise_permutations[, 1],
assoc_protein = pairwise_permutations[, 2]
) %>%
mutate(
shared_pathway = map2(`protein`, `assoc_protein`, function(x, y) {
per_permutation(pathway_df = pathway_df, permutation = c(x, y))
}),
combination = map2_chr(`protein`, `assoc_protein`, function(x, y) {
comb <- sort(c(x, y))
paste(comb, collapse = ",")
}),
count = map_int(`shared_pathway`, function(x) na.omit(x) %>% length()),
connected = count > 0
)
shared_pathway_df
}
shared_pathway_df <- get_shared_pathway(pathway_df)
n_pairs <- length(shared_pathway_df %>% filter(count > 0))
print(glue::glue("Num. shared_pathway pairs: {n_pairs}"))
#> Num. shared_pathway pairs: 6
shared_pathway_df %>% arrange(desc(count))
#> # A tibble: 72 x 6
#> protein assoc_protein shared_pathway combination count connected
#> <chr> <chr> <list> <chr> <int> <lgl>
#> 1 P05362 Q9UMF0 <chr [2]> P05362,Q9UMF0 2 TRUE
#> 2 Q9UMF0 P05362 <chr [2]> P05362,Q9UMF0 2 TRUE
#> 3 P05362 P29597 <chr [1]> P05362,P29597 1 TRUE
#> 4 P29597 P05362 <chr [1]> P05362,P29597 1 TRUE
#> 5 O75791 P05362 <chr [0]> O75791,P05362 0 FALSE
#> 6 O75791 P0C6A0 <NULL> O75791,P0C6A0 0 FALSE
#> 7 O75791 P29597 <chr [0]> O75791,P29597 0 FALSE
#> 8 O75791 Q13445 <NULL> O75791,Q13445 0 FALSE
#> 9 O75791 Q6P4F2 <chr [0]> O75791,Q6P4F2 0 FALSE
#> 10 O75791 Q8N9T8 <NULL> O75791,Q8N9T8 0 FALSE
#> # … with 62 more rowsWe can further query EpiGraphDB regarding the detailed pathway information using GET /meta/nodes/Pathway/search.
get_pathway_info <- function(reactome_id) {
endpoint <- "/meta/nodes/Pathway/search"
params <- list(id = reactome_id)
df <- query_epigraphdb(route = endpoint, params = params, mode = "table")
df
}
pathway <- shared_pathway_df %>%
pull(shared_pathway) %>%
unlist() %>%
unique()
pathway_info <- pathway %>% map_df(get_pathway_info)
pathway_info %>% print()
#> # A tibble: 3 x 3
#> node.in_disease node.name node.reactome_id
#> <chr> <chr> <chr>
#> 1 FALSE Interleukin-4 and Interleukin-13 signaling R-HSA-6785807
#> 2 FALSE Integrin cell surface interactions R-HSA-216083
#> 3 FALSE Immunoregulatory interactions between a Lymp… R-HSA-198933In order to extract protein groups from the shared pathways, the last step for this query is to convert the shared pathway data into a graph where
- Proteins are represented as nodes.
- And shared pathways are represented as edges.
- Then the protein groups (via shared pathways) are the connected communities in the graph.
We then count the number of nodes in each connected community and plot the graph.
protein_df_to_graph <- function(df) {
df_connected <- df %>%
filter(connected) %>%
distinct(`combination`, .keep_all = TRUE)
nodes <- df %>%
pull(protein) %>%
unique()
graph <- igraph::graph_from_data_frame(
df_connected,
directed = FALSE, vertices = nodes
)
graph$layout <- igraph::layout_with_kk
graph
}
graph_to_protein_groups <- function(graph) {
graph %>%
igraph::components() %>%
igraph::groups() %>%
tibble(group_member = .) %>%
mutate(group_size = map_int(`group_member`, length)) %>%
arrange(desc(group_size))
}
pathway_protein_graph <- shared_pathway_df %>% protein_df_to_graph()
pathway_protein_groups <- pathway_protein_graph %>% graph_to_protein_groups()
pathway_protein_groups %>% str()
#> tibble [7 × 2] (S3: tbl_df/tbl/data.frame)
#> $ group_member:List of 7
#> ..$ 2: chr [1:3] "P05362" "P29597" "Q9UMF0"
#> ..$ 1: chr "O75791"
#> ..$ 3: chr "P0C6A0"
#> ..$ 4: chr "Q13445"
#> ..$ 5: chr "Q6P4F2"
#> ..$ 6: chr "Q8N9T8"
#> ..$ 7: chr "Q9BYD3"
#> ..- attr(*, "dim")= int 7
#> ..- attr(*, "dimnames")=List of 1
#> .. ..$ : chr [1:7] "2" "1" "3" "4" ...
#> $ group_size : Named int [1:7] 3 1 1 1 1 1 1
#> ..- attr(*, "names")= chr [1:7] "2" "1" "3" "4" ...plot(pathway_protein_graph)
