本文是洲更寫的，文中說到兩步，其實只需要這一步，然后就可以用enricher或GSEA函數(shù)進行分析，也就是通用的分析函數(shù)，不管是什么注釋，不管是什么物種，都可以支持。 這樣做的好處也有一些 第二步所謂需要序列，其實是想本地blast，自己注釋，這種一般是和KO數(shù)據(jù)庫進行比對。等候洲更的更新，看他后面怎么玩。
Y叔的clusterProfiler的一大優(yōu)點就是能夠利用最新的KEGG數(shù)據(jù)庫，而不是停留在最后一個公開版的KEGG數(shù)據(jù)庫(2011-5-15). 大部分情況下，大家都是直接用enrichKEGG()或者gseKEGG()完成富集分析，但是我最近想到，我其實可以利用這個功能反向建立一個注釋用的數(shù)據(jù)庫，只需要兩步 第一步: 獲取KEGG編號對應的基因編號
第二步: 根據(jù)基因編號獲取序列信息 這里只講第一步，如何用Y叔clusterProfiler獲取最新的KEGG和基因的對應關(guān)系， 我們以人類為例。
先用download_KEGG下載給定物種的KEGG數(shù)據(jù)庫 hsa_kegg <- clusterProfiler::download_KEGG('hsa')
這一步得到的是一個列表，這個列表有兩個成員， names(hsa_kegg)
[1] 'KEGGPATHID2EXTID' 'KEGGPATHID2NAME'
一個是KEGG的通路編號和基因編號的關(guān)系，另一個是KEGG通路編號和名字的關(guā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的編號？這可以使用biomaRt進行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的基因編號，通過merge就可以將其添加到之前的數(shù)據(jù)框中 PATH_ID_NAME <- merge(PATH_ID_NAME, ensembl_gene_id, by.x= 'ENTREZID',by.y= 'entrezgene')
|
