使用clusterProfiler進行KEGG富集分析

生信修煉手冊 >《待分類》

2019.12.24

關(guān)注

KEGG pathway是最常用的功能注釋數(shù)據(jù)庫之一，可以利用KEGG 的API獲取一個物種所有基因?qū)?yīng)的pathway注釋，human對應(yīng)的API 鏈接如下

http://rest.kegg.jp/link/hsa/pathway

通過該鏈接可以獲得以下內(nèi)容

path:hsa00010 hsa:10327path:hsa00010 hsa:124path:hsa00010 hsa:125

第一列為pathway編號，第二列為基因編號。這里只提供了pathway編號，我們還需要pathway對應(yīng)的描述信息，同樣也可以通過以下API鏈接得到

http://rest.kegg.jp/list/

通過該鏈接可以獲得如下內(nèi)容

path:map00010 Glycolysis / Gluconeogenesispath:map00020 Citrate cycle (TCA cycle)path:map00030 Pentose phosphate pathwaypath:map00040 Pentose and glucuronate interconversionspath:map00051 Fructose and mannose metabolism

第一列為pathway編號，第二列為具體的描述信息。需要注意的是，pathway是一個跨物種的概念，原始的pathway編號為map或者ko加數(shù)字，對于特定物種，改成物種對應(yīng)的三字母縮寫, 比如human對應(yīng)hsa, 所有擁有pathway信息的物種和對應(yīng)的三字母縮寫見如下鏈接

https://www.genome.jp/kegg/catalog/org_list.html

clusterProfiler也是通過KEGG API去獲取物種對應(yīng)的pathway注釋，對于已有pathway注釋的物種，我們只需要知道對應(yīng)的三字母縮寫， clusterProfiler就會聯(lián)網(wǎng)自動獲取該物種的pathway注釋信息。

和GO富集分析類似，對于KEGG的富集分析也包含以下兩種

1. Over-Representation Analysis

過表達分析其實就是費舍爾精確檢驗，分析的代碼如下

ego <- enrichKEGG(  gene          = gene,  keyType     = "kegg",  organism   = 'hsa',  pvalueCutoff      = 0.05,  pAdjustMethod     = "BH",  qvalueCutoff  = 0.05)

我們只需要提供差異基因的列表和物種對應(yīng)的三字母縮寫，默認基因ID為kegg gene id, 通過keyType參數(shù)指定，也可以是ncbi-geneid, ncbi-proteind, uniprot。

不同類型ID的轉(zhuǎn)換也是通過KEGG API實現(xiàn)的，比如hsa的kegg gene id和ncbi-geneid的對應(yīng)關(guān)系見以下鏈接

http://rest.kegg.jp/conv/ncbi-geneid/hsa

hsa:1 ncbi-geneid:1hsa:100009667 ncbi-geneid:100009667hsa:100009676 ncbi-geneid:100009676hsa:10 ncbi-geneid:10hsa:100 ncbi-geneid:100

在clusterProfiler中，可以通過bitr_kegg函數(shù)，調(diào)用KEGG API, 來實現(xiàn)ID 轉(zhuǎn)換功能，示例如下

bitr_kegg("1",fromType = "kegg",toType = 'ncbi-proteinid',organism='hsa')

2. Gene Set Enrichment Analysis

對應(yīng)的函數(shù)為gseKEGG, 用法如下

kk <- gseKEGG(  geneList  = gene,  keyType  = 'kegg',  organism = 'hsa',  nPerm  = 1000,  minGSSize = 10,  maxGSSize = 500,  pvalueCutoff = 0.05,  pAdjustMethod     = "BH")

對于pathway數(shù)據(jù)庫中沒有的物種，也支持讀取基因的pathway注釋文件，然后進行分析，注釋文件的格式如下

GeneId	KEGG	Description
1	ko:00001	spindle
2	ko:00002	mitotic spindle
3	ko:00003	kinetochore

只需要3列信息即可，第一列為geneID, 第二列為基因?qū)?yīng)的pathway編號，第三列為pathway的描述信息。這3列的順序是無所謂的，只要包含這3種信息就可以了。

讀取該文件，進行分析的代碼如下

data <- read.table(  "pathway_annotation.txt",  header = T,  sep = "\t")go2gene <- data[, c(2, 1)]go2name <- data[, c(2, 3)]# 費舍爾精確檢驗x <- enricher(gene,TERM2GENE = go2gene,TERM2NAME = go2name)# GSEA富集分析x <- GSEA(gene,TERM2GENE = go2gene,TERM2NAME = go2name)

對于KEGG富集分析的結(jié)果，clusterProfiler提供了以下幾種可視化策略

1. barplot

用散點圖展示富集到的pathways,用法如下

barplot(kk, showCategory = 10)

生成的圖片如下

橫軸為該pathway的差異基因個數(shù)，縱軸為富集到的pathway的描述信息， showCategory指定展示的pathway的個數(shù)，默認展示顯著富集的top10個，即p.adjust最小的10個。注意的顏色對應(yīng)p.adjust值，從小到大，對應(yīng)藍色到紅色。