Y叔的clusterProfiler的一大优点就是能够利用最新的KEGG数据库，而不是停留在最后一个公开版的KEGG数据库(2011-5-15).
大部分情况下，大家都是直接用enrichKEGG()或者gseKEGG()完成富集分析，但是我最近想到，我其实可以利用这个功能反向建立一个注释用的数据库，只需要两步
第一步: 获取KEGG编号对应的基因编号
第二步: 根据基因编号获取序列信息
这里只讲第一步，如何用Y数clusterProfiler获取最新的KEGG和基因的对应关系， 我们以人类为例。
先用download_KEGG下载给定物种的KEGG数据库
hsa_kegg <- clusterProfiler::download_KEGG("hsa")
这一步得到的是一个列表，这个列表有两个成员，
names(hsa_kegg)
[1] "KEGGPATHID2EXTID" "KEGGPATHID2NAME"
一个是KEGG的通路编号和基因编号的关系，另一个是KEGG通路编号和名字的关系
我们可以将其进行合并
PATH2ID <- hsa_kegg$KEGGPATHID2EXTID
PATH2NAME <- hsa_kegg$KEGGPATHID2NAME
PATH_ID_NAME <- merge(PATH2ID, PATH2NAME, by="from")
colnames(PATH_ID_NAME) <- c("KEGGID", "ENTREZID", "DESCRPTION")
最后的PATH_ID_NAME的表格信息如下
KEGG关系表
保存到本地
write.table(PATH_ID_NAME, "HSA_KEGG.txt", sep="\t")
拓展：如何在此基础上增加ENSEMBL的编号？这可以使用biomaRt进行ID转换
library(biomaRt)
mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl"))
entrezgene <- PATH_ID_NAME$ENTREZID
# This step need some time
ensembl_gene_id<- getBM(attributes=c("ensembl_gene_id", "entrezgene"),
filters = "entrezgene",
values=entrezgene , mart= mart)
于是我们得到了ENSEMBL的基因编号，通过merge就可以将其添加到之前的数据框中
PATH_ID_NAME <- merge(PATH_ID_NAME, ensembl_gene_id, by.x= "ENTREZID",by.y= "entrezgene")
本文还可在http://xuzhougeng.top/阅读
版权声明：本博客所有文章除特别声明外，均采用 知识共享署名-非商业性使用-禁止演绎 4.0 国际许可协议 (CC BY-NC-ND 4.0) 进行许可。
扫码即刻交流
