本文是洲更寫的，文中說(shuō)到兩步，其實(shí)只需要這一步，然后就可以用
enricher或GSEA函數(shù)進(jìn)行分析，也就是通用的分析函數(shù)，不管是什么注釋，不管是什么物種，都可以支持。這樣做的好處也有一些
分析時(shí)不用聯(lián)網(wǎng)，可以搞個(gè)local版本
分析具有可重復(fù)性，實(shí)時(shí)爬網(wǎng)絡(luò)數(shù)據(jù)，畢竟有時(shí)候服務(wù)器端有更新的話，就會(huì)導(dǎo)致有一些變化。
第二步所謂需要序列，其實(shí)是想本地blast，自己注釋，這種一般是和KO數(shù)據(jù)庫(kù)進(jìn)行比對(duì)。等候洲更的更新，看他后面怎么玩。
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叔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)，通過(guò)merge就可以將其添加到之前的數(shù)據(jù)框中
PATH_ID_NAME <- merge(PATH_ID_NAME, ensembl_gene_id, by.x= 'ENTREZID',by.y= 'entrezgene')
聯(lián)系客服
