r 从KEGG信号传导途径中拉出所有磷酸化边缘

Posted 2021-05-24

require(KEGGgraph, quietly = TRUE)
require(dplyr, quietly = TRUE)
require(magrittr, quietly = TRUE)
# install.packages("devtools")
#devtools::install_github("robertness/lucy")
require(lucy)
maps <- c(mapk = "04010", pi3 = "04151", tcell = "04660", ras = "04014", rap1 = "04015",
          erb = "04012", wnt = "04310", tgfb = "04350", hippo = "04390", vegf = "04370", 
          jakstat = "04630", nfkb = "04064", tnf = "04668", hif1 = "04066", foxo = "04068",
          calcium = "04020", camp = "04024", cgmp = "04022", ampk = "04152", mtor = "04150")
vertex_master <- NULL
edge_master <- NULL
for(map in maps){
  g_nell <- tempfile() %T>%
    {retrieveKGML(map, organism="hsa", destfile=., method="curl", quiet=TRUE)} %>%
    parseKGML2Graph(expandGenes=FALSE) 
  vertex_list <- getKEGGnodeData(g_nell) %>%
    {data.frame(
      kegg = unlist(lapply(., function(item) item@name[1])),
      label = unlist(lapply(., function(item)
        strsplit(item@graphics@name, ",")[[1]][1])), stringsAsFactors = F)}
  g_init <- igraph.from.graphNEL(g_nell) 
  V(g_init)$name <- vertex_list$kegg 
  vertex_list <- filter(vertex_list, !duplicated(kegg))
  edge_list <- getKEGGedgeData(g_nell) %>%
    lapply(function(item){
      if(length(item@subtype) > 0){
        subtype_info <- item@subtype
        # KEGG uses a hierarchy of term for describing terms
        # for example, the first edge type is "activation", the second is "phosphorylation"
        # where phosphorylation is a type of activation.  The second term is more specific than
        # the first, so when it is provided, use it in lieu of the first type.
        if(length(subtype_info) > 1) {
          return(subtype_info[[2]]@name)
        } else {
          return(subtype_info$subtype@name)
        }
      } 
      NA
    }) %>%
    unlist %>%
    {cbind(get.edgelist(g_init), type = .)} %>%
    data.frame %>%
    filter(type == "phosphorylation")
  edge_master <- rbind(edge_master, edge_list)
  vertex_master <- rbind(vertex_master, vertex_list)
}
edge_master <- edge_master %>%
  as.data.frame %>%
  unique
vertex_master <- vertex_master %>%
  unique %>%
  filter(!duplicated(kegg))
g <- graph.data.frame(edge_master, directed = TRUE, vertices = vertex_master)
V(g)$kid <- V(g)$name
V(g)$name <- V(g)$label
g <- g - V(g)[igraph::degree(g) == 0]