Thank you for your answer, I managed to adapt it to Pathways instead of Modules :
library(tidyverse)
library(KEGGREST)
kegg_id_names = c('hsa:7040')
path_ids = map(kegg_id_names, ~keggLink('pathway', .)) %>%
unlist() %>%
unique()
path_data = keggGet(path_ids[1:3])
get_kegg_position = function(path_record) {
gene_ids = path_record$GENE[seq(1,length(path_record$GENE),2)]
gene_names = gsub(";.*","",path_record$GENE[seq(2,length(path_record$GENE),2)])
if (is.null(gene_ids)) {
res = tibble(gene = NA,
position = NA)
}else{
gene_posns = seq(1, 0, length.out = length(gene_ids))
res = tibble(gene_id = gene_ids,
gene_name = gene_names,
position = gene_posns) %>%
separate_rows('gene_id', sep = ',')
}
res %>% mutate(path_id = path_record$ENTRY,
path_name = path_record$NAME)
}
kegg_positions_by_module = map_df(path_data,
get_kegg_position) %>%
filter(!is.na(position))
# TGFB1 score in MAPK pathway (should be close to 1)
kegg_positions_by_module %>%
filter(path_id=="hsa04010") %>% filter(gene_name == "TGFB1")
# SRF score in MAPK pathway (should be lower than TGFB1 as it is only upstream of 1 gene (FOS))
kegg_positions_by_module %>%
filter(path_id=="hsa04010") %>% filter(gene_name == "SRF")
# MAX score in MAPK pathway (should be close to 0)
kegg_positions_by_module %>%
filter(path_id=="hsa04010") %>% filter(gene_name == "MAX")
Here I obtain a value of only 0.433 for TGFB1 whereas I should get value closer to 1. For SRF gene, I obtain value of 0.481 which is higher than TGFB1 while SRF is upstream of only 1 gene and TGFB1 is upstream of many more genes. (See MAPK pathway diagram)
I think your approach is based on the hypothesis that KEGG ranks genes from more upstream to more downstream. While this might be true for reactions, I think it is not for Pathways. I think the order of genes in KEGG pathway is from Left to Right and Top to Bottom and in the KEGG diagrams.
I am sure this score would never be accurately quantitative but I would hope it gives an idea about how often a gene is upstream of pathways.
Hello Alex,
I used KEGGgraph as suggested to get Directed Graphs from KEGG pathway database for each gene in my gene list. Here is a solution, the code is probably suboptimal but I hope it is understandable.
The first fonction
GetKEGGigraphtakes in a KEGG pathway id (e.g. "04010" for MAPK signalling pathway). It first download the pathway as Directed Graph using KEGGgraph, then transform it into an "igraph graph" object. Each gene is a node and each relation geneA -> geneB is a directed edge. We first remove isolated nodes & self loops from the graphs. Then, as some pathways are Cyclic, we need to make them acyclic using Minimum Spanning Tree algorithm. This is because we cannot make a proper gene hierarchy in cyclic graphs. TheGetKEGGigraphfunction returns the Minimum Spanning Tree.The second function
upstream_score_KEGGtakes as input a gene list as gene names (e.g.c("TGFB1","MAX")). It first retrieve for each gene all the pathways containing the gene. Then, for each gene in each pathway, it callsGetKEGGigraphto get the Minimum Spanning Tree. The tree is sorted using 'Topological Sorting' algorithm (From wikipedia :"In computer science, a topological sort or topological ordering of a directed graph is a linear ordering of its vertices such that for every directed edge uv from vertex u to vertex v, u comes before v in the ordering"). This way, most upstream genes will be placed before more downstream genes. The functionupstream_score_KEGGoutputs a data.frame with the upstream score (1 - relative position) of the gene in each of it's pathway, or NA if no pathway is found.I have very little experience with graph theory so if I would gladly take advices on this part, particularly : it this the best way to get an 'order' of genes in a graph ? I think that with topological sorting, two 'unrelated nodes' will be placed in a kind of arbitrary way.
Here is the code:
So the results is, for MAPK signalling pathway, TGFB1 has score of 0.731 & MAX has score of 0.187. This gives you an idea of how upstream the gene is in a pathway. In the 'TGF-beta signaling pathway' for instance the score for TGFB1 is only 0.31 because it is downstream of 5 genes in this pathway (TGF-beta signaling pathway).
The trouble comes when you average over all the pathways, because as you said pathways are complex and the end of any pathways usually comes before and after any pathway.... So the same gene might be the most downstream in one pathway and the most upstream in another pathway...
I hope the approach is not too complex for nothing and I would gladly hear any comments on better way to achieve the same results.
I must put the functions separately as the limit is 5,000 characers. Here is the first function
Here is the second function:
You can put all this code in a GitHub
gist. Make it public and then paste that link into your comment above. Biostars code will then render it all in one stretch.Please use
ADD COMMENT/ADD REPLYwhen responding to existing posts to keep threads logically organized.SUBMIT ANSWERis for new answers to original question.Okay I will do that, thank you
If you can delete these duplicate comments that would be great. You can find "delete post" option under the
moderatebutton for your own posts.