我们经常要用到KEGG数据库来对基因做功能分析。经常有长得好看的朋友问：如何获得整个通路的基因信息？
其实我们有多种方法可以获得通路中所有的基因情况，本文通过KEGG的原始网页生成某个通路的基因表格。
准备文件：
image.png
到KEGG网站获取某个通路的KGML文件（Download KGML），下载下来的其实是一个xml文件。看一下xml的结构：
image.png
我们想要的信息都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)
接着就用我们之前文章：R语者小case之——从GTF文件生成注释表格做基因ID转换 所获得的注释表格，把这里的name转换为其他gene ID，生成一个完整的表格。
