R语者小case之——从KEGG原始网页批量生成通路的基因表格

本文链接：https://blog.csdn.net/weixin_49352797/article/details/118672140

我们经常要用到KEGG数据库来对基因做功能分析。经常长得好看的朋友问：如何获得整个通路的基因？

其实我们有多种方法可以获得通路中所有的基因情况，本文通过KEGG的原始网页生成某个通路的基因表格。

准备文件：
在这里插入图片描述
到KEGG网站获取某个通路的KGML文件（Download KGML），下载下来的其实是一个xml文件。看一下xml的结构：

我们想要的信息都entry这一行，之后的事情就交给R了。

library(plyr)#需要载入这两个包，没有的话，请先安装
library(XML)
readin <- xmlToList("hsa04060.xml")#将xml转为list
readin1 <- readin[grepl("^entry",names(readin))]#保留entry开头的行
#readin2 <- readin[!grepl("^entry",names(readin))]
sub_id <- lapply(readin1,function(x){
  x2 <- x$.attrs
  table1 <- data.frame(t(x2))
  rownames(table1) <- table1$id
  return(table1)
})
#用上了lapply函数，批量运行，将entry行的内容保存为data.frame
sub_tb <- do.call(rbind.fill,sub_id)
#再用do.call做了一个批量rbind，这里用rbind.fill自动补充缺值为NA
head(sub_tb)
  id      name type                                             link
1 50 hsa:53833 gene https://www.kegg.jp/dbget-bin/www_bget?hsa:53833
2 51 hsa:53833 gene https://www.kegg.jp/dbget-bin/www_bget?hsa:53833
3 52  hsa:3565 gene  https://www.kegg.jp/dbget-bin/www_bget?hsa:3565
4 53   hsa:659 gene   https://www.kegg.jp/dbget-bin/www_bget?hsa:659
5 54    hsa:93 gene    https://www.kegg.jp/dbget-bin/www_bget?hsa:93
6 55    hsa:91 gene    https://www.kegg.jp/dbget-bin/www_bget?hsa:91
tail(sub_tb)#看到缺值被自动填充为NA了
     id      name  type link
436 718 undefined group <NA>
437 719 undefined group <NA>
438 720 undefined group <NA>
439 721 undefined group <NA>
440 722 undefined group <NA>
441 723 undefined group <NA>

这里没有用循环语句来实现，而用lapply，do.call做批量运行，用好这一系列函数，事半功倍。

这个时候，有些颜值高的朋友可能会问：这只是其中一个通路的，我要所有通路的该怎么做？

安排！在KEGG的API安排一个“小爬虫”。

  library(curl)#需要载入这两个包，没有的话，请先安装
  library(stringr)
 pathway_all <- read.table("http://rest.kegg.jp/list/pathway/hsa",sep = "\t")
  #从KEGG获取人所有的通路信息
  pathID <- str_split_fixed(pathway_all$V1,":",2)[,2]
  #从上表获得所有pathwayID
  for(i in 1:length(pathID)){
    pathurl <- paste0("http://rest.kegg.jp/get/",pathID[i],"/kgml")
    dir.create("./xml/")
    xml_file <- paste0("./xml/",pathID[i],".xml")
    curl_download(pathurl, xml_file)#用小爬虫把数据爬下来
  }
  library(plyr)#需要载入这两个包，没有的话，请先安装
  library(XML)
  sub_ls <- list()
  for(pathID in pathID){
    readin <- xmlToList(paste0("./xml/",pathID,".xml"))#将xml转为list
    readin1 <- readin[grepl("^entry",names(readin))]#保留entry开头的行
    #readin2 <- readin[!grepl("^entry",names(readin))]
    sub_id <- lapply(readin1,function(x){
      x2 <- x$.attrs
      table1 <- data.frame(t(x2))
      rownames(table1) <- table1$id
      return(table1)
    })
    #用上了lapply函数，批量运行，将entry行的内容保存为data.frame
    sub_ls[[pathID]] <- do.call(rbind.fill,sub_id)
    #再用do.call做了一个批量rbind，这里用rbind.fill自动补充缺值为NA
  }
  sub_tb <- do.call(rbind.fill,sub_ls)