使用R語言clusterProfiler包進(jìn)行KEGG富集分析（有參分析）
我們在前面的視頻教程中取胎，為大家講述了GO展哭、KEGG富集分析的基本原理湃窍，并展示了如何使用DAVID進(jìn)行GO富集在線分析的過程，詳情點擊查看“沒想到基因GO富集分析這么簡單”匪傍。本節(jié)繼續(xù)以R語言clusterProfiler包中的方法為例您市，展示如何進(jìn)行KEGG富集分析。
圖文教程
以下是圖文版的簡化教程役衡。
1 準(zhǔn)備輸入數(shù)據(jù)
我們提供的示例數(shù)據(jù)為人類基因茵休，文件中共包含兩列。第一列是基因的ENTREZ ID手蝎，用于和KEGG數(shù)據(jù)庫建立聯(lián)系對應(yīng)基因以便進(jìn)行富集分析榕莺；第二列是基因的SYMBOL ID，也就是通俗的基因名稱棵介，用于基因名稱轉(zhuǎn)換用钉鸯，詳見下文過程或上文視頻。
注：clusterProfiler默認(rèn)只能識別3種類型的基因名稱用作富集分析邮辽，ENTREZ ID唠雕、NCBI ID以及UniProt ID。我們的示例文件中就使用的ENTREZ ID吨述。如果老師或同學(xué)們不清楚怎樣轉(zhuǎn)換基因ID岩睁，例如ENSEMBL ID向ENTREZ ID的轉(zhuǎn)換，可以添加微信公眾號揣云，聯(lián)系我們尋求幫助捕儒。
2 clusterProfiler的KEGG富集分析
讀取基因列表文件至R中，加載并使用clusterProfiler的內(nèi)部函數(shù)enrichKEGG()即可自動完成KEGG富集分析邓夕。
注：以下方法均對于KEGG數(shù)據(jù)庫中已經(jīng)收錄的物種而言的刘莹，若目標(biāo)物種未包含在KEGG數(shù)據(jù)庫中，則無法進(jìn)行KEGG富集翎迁。例如該示例的物種是人類栋猖，其它物種查詢可點擊以下鏈接：
https://www.kegg.jp/kegg-bin/get_htext?htext=br08601_ko00001&query=ko&option=-w
#bioconductor 安裝 clusterProfiler
#install.packages('BiocManager') #需要首先安裝 BiocManager，如果尚未安裝請先執(zhí)行該步
#BiocManager::install('clusterProfiler')
?
#加載R包
library(clusterProfiler)
?
#讀取待富集的基因名稱（示例數(shù)據(jù)汪榔，基因名稱列表）
gene_list <- read.delim('gene_list.txt', stringsAsFactors = FALSE)
?
#每次打開R計算時蒲拉，他會自動連接kegg官網(wǎng)獲得最近的物種注釋信息肃拜，因此數(shù)據(jù)庫一定都是最新的
kegg <- enrichKEGG(
gene = gene_list$ENTREZ, #指定基因名稱，例如 ENTREZ ID
keyType = 'kegg', #指定基因名稱類型雌团，kegg 代表了使用 ENTREZ ID
organism = 'hsa', #hsa代表人類燃领，其它物種更改這行即可椿猎，可在 KEGG 官網(wǎng)查詢物種名縮寫
pvalueCutoff = 1, #p 值 p 調(diào)整值之類的脖律，指定 1 為閾值，也就是將所有的條目都輸出了
qvalueCutoff = 1,
pAdjustMethod = 'fdr' #p 值調(diào)整方法
)
?
#輸出旗扑，默認(rèn)情況下基因名稱顯示為富集時使用的名稱
write.table(kegg, 'KEGG.ENTREZ.txt', sep = '\t', quote = FALSE, row.names = FALSE)
最后輸出富集分析結(jié)果灵寺，對于各列內(nèi)容：
ID和Description曼库，富集到的KEGG條目的名稱和描述；
GeneRatio和BgRatio略板，分別為富集到該KEGG條目中的基因數(shù)目/給定基因的總數(shù)目毁枯，以及該條目中背景基因總數(shù)目/該物種所有已知的KEGG功能基因數(shù)目；
pvalue叮称、p.adjust和qvalue种玛，富集p值、校正后p值和q值信息瓤檐；
geneID和Count赂韵，分別為富集到該KEGG條目中的基因名稱和數(shù)目，由于分析中使用的ENTREZ ID挠蛉，故這里也默認(rèn)顯示為ENTREZ ID祭示。
3 基因名稱轉(zhuǎn)換
由于ENTREZ ID只是一串?dāng)?shù)字編號，不利于辨認(rèn)碌秸，因此我們可以考慮將它轉(zhuǎn)換為通俗的基因名稱绍移，也就是使用SYMBOL ID。
#基因名稱轉(zhuǎn)換讥电，例如將基因 ENTREZ ID 轉(zhuǎn)換為 SYMBOL ID
kegg <- read.delim('KEGG.ENTREZ.txt', stringsAsFactors = FALSE)
rownames(gene_list) <- as.character(gene_list$ENTREZ)
?
for (i in 1:nrow(kegg)) {
new <- c()
old <- unlist(strsplit(kegg[i,'geneID'], '/'))
for (id in old) new <- c(new, gene_list[id,'SYMBOL'])
kegg[i,'geneID'] <- paste(new, collapse = '/')
}
?
write.table(kegg, 'KEGG.SYMBOL.txt', sep = '\t', quote = FALSE, row.names = FALSE)
新輸出的表格中蹂窖，geneID中的基因名稱就轉(zhuǎn)換為SYMBOL ID了，這樣我們就可以很清楚地知道富集條目中所包含的基因具體是哪種類型的恩敌。
4 clusterProfiler附帶的作圖功能
此外瞬测，clusterProfiler中也額外提供了一系列的可視化方案用于展示本次富集分析結(jié)果，具有極大的便利纠炮。
#clusterProfiler 包里的一些默認(rèn)作圖方法月趟，例如
barplot(kegg) #富集柱形圖
dotplot(kegg) #富集氣泡圖
cnetplot(kegg) #網(wǎng)絡(luò)圖展示富集功能和基因的包含關(guān)系
emapplot(kegg) #網(wǎng)絡(luò)圖展示各富集功能之間共有基因關(guān)系
heatplot(kegg) #熱圖展示富集功能和基因的包含關(guān)系
這樣，KEGG富集分析結(jié)果表連同富集柱形圖恢口、氣泡圖等一起就都獲得了孝宗，是不是很快捷方便呢？
