Y叔的clusterProfiler的一大優(yōu)點(diǎn)就是能夠利用最新的KEGG數(shù)據(jù)庫(kù)怔檩,而不是停留在最后一個(gè)公開版的KEGG數(shù)據(jù)庫(kù)(2011-5-15).
大部分情況下磷瘤,大家都是直接用enrichKEGG()
或者gseKEGG()
完成富集分析,但是我最近想到,我其實(shí)可以利用這個(gè)功能反向建立一個(gè)注釋用的數(shù)據(jù)庫(kù),只需要兩步
第一步: 獲取KEGG編號(hào)對(duì)應(yīng)的基因編號(hào)
第二步: 根據(jù)基因編號(hào)獲取序列信息
這里只講第一步认烁,如何用Y數(shù)clusterProfiler獲取最新的KEGG和基因的對(duì)應(yīng)關(guān)系, 我們以人類為例介汹。
先用download_KEGG
下載給定物種的KEGG數(shù)據(jù)庫(kù)
hsa_kegg <- clusterProfiler::download_KEGG("hsa")
這一步得到的是一個(gè)列表却嗡,這個(gè)列表有兩個(gè)成員,
names(hsa_kegg)
[1] "KEGGPATHID2EXTID" "KEGGPATHID2NAME"
一個(gè)是KEGG的通路編號(hào)和基因編號(hào)的關(guān)系嘹承,另一個(gè)是KEGG通路編號(hào)和名字的關(guān)系
我們可以將其進(jìn)行合并
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的表格信息如下
保存到本地
write.table(PATH_ID_NAME, "HSA_KEGG.txt", sep="\t")
拓展:如何在此基礎(chǔ)上增加ENSEMBL的編號(hào)窗价?這可以使用biomaRt進(jìn)行ID轉(zhuǎn)換
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的基因編號(hào),通過merge就可以將其添加到之前的數(shù)據(jù)框中
PATH_ID_NAME <- merge(PATH_ID_NAME, ensembl_gene_id, by.x= "ENTREZID",by.y= "entrezgene")
本文還可在http://xuzhougeng.top/閱讀
版權(quán)聲明:本博客所有文章除特別聲明外叹卷,均采用 知識(shí)共享署名-非商業(yè)性使用-禁止演繹 4.0 國(guó)際許可協(xié)議 (CC BY-NC-ND 4.0) 進(jìn)行許可撼港。