尋找差異表達(dá)的基因并識(shí)別它們的功能掷豺，是我們進(jìn)行RNA測(cè)序的最主要目的默辨。很明顯表谊，這些差異的基因必然與功能改變密切相關(guān)，例如，比較患病個(gè)體與正常個(gè)體的組織表達(dá)譜，不難想到這些顯著失調(diào)的基因參與了生物學(xué)過程雾消、信號(hào)通路等，導(dǎo)致了疾病的發(fā)生破讨。
前面已經(jīng)講了如何使用DESeq2隙笆、edgeR基于轉(zhuǎn)錄組測(cè)序獲得的基因表達(dá)值鑒定差異表達(dá)基因仰冠。那么，后續(xù)如何繼續(xù)通過生信分析的方法，探索差異表達(dá)的基因發(fā)揮了怎樣的功能，參與了哪些調(diào)控通路呢银择？
我們平時(shí)看RNA-seq相關(guān)的文獻(xiàn)時(shí)，文章中在鑒定了差異表達(dá)的基因后袜瞬，大都會(huì)在隨后承接幾句關(guān)于這些失調(diào)基因所涉及通路的描述。例如，討論這些差異基因主要映射到哪些GO或KEGG分類條目中，以說明基因表達(dá)的改變會(huì)導(dǎo)致哪些調(diào)控途徑原有功能失調(diào)逐哈，進(jìn)而與表型聯(lián)系起來。通常稱這種分析為GO蚀腿、KEGG富集分析。
目前瞎嬉，能夠進(jìn)行GO便监、KEGG富集分析的工具有很多逊移，不同的工具之間在算法、數(shù)據(jù)庫組成上略有不同案铺，因此結(jié)果也可能大相徑庭暇番。本篇就先以R包c(diǎn)lusterProfiler的方法為例岳服，展示如何基于給定的基因列表分析它們的GO璃搜、KEGG功能唾糯。
clusterProfiler包的安裝
對(duì)于clusterProfiler的安裝也很簡(jiǎn)單，一般情況下脐帝，直接通過Bioconductor安裝clusterProfiler就可以了禁偎。
#bioconductor 安裝
#install.packages('BiocManager') #需要首先安裝 BiocManager，如果尚未安裝請(qǐng)先執(zhí)行該步
BiocManager::install('clusterProfiler')
clusterProfiler的GO富集分析（有參向）
首先來看GO富集分析。
注：這里均對(duì)于有參考基因組的情況而言的，無參分析暫不涉及。
1牲蜀、準(zhǔn)備輸入數(shù)據(jù)
待分析的數(shù)據(jù)就是一串基因名稱了，可以是ensembl id、entrze id或者symbol id等類型都可以先慷。把基因名稱以一列的形式排開无午，放在一個(gè)文本文件中（例如命名“gene.txt”）踩验。Excel中查看钠龙，就是如下示例這種樣式咬腋。
2撒汉、加載參考物種的基因注釋數(shù)據(jù)庫
對(duì)于有參考基因組物種的分析，首先需要指定該物種的基因數(shù)據(jù)庫。
存在兩組情況寻歧，一種是常見物種噪珊，另一種是不常見物種阵难。
#（1）對(duì)于常見物種播玖，例如人類果覆，有些專門的 R 包數(shù)據(jù)庫，例如人類參考基因組 hg19 的
library(org.Hs.eg.db)
#（2）對(duì)于不常見的物種，但卻是存在參考基因組的情況
#通過 AnnotationHub 包索引基因組，例如期望找綿羊（Ovis aries）的注釋庫
library(AnnotationHub)
hub <- AnnotationHub()
query(hub, 'Ovis aries') #輸入綿羊（Ovis aries）的名稱進(jìn)行匹配
sheep <- hub[['AH72269']] #返回了數(shù)據(jù)庫編號(hào) AH72269倒淫，就可以加載該庫
3兴枯、GO富集分析
加載了注釋庫之后躺坟，讀取基因列表文件，并使用clusterProfiler的內(nèi)部函數(shù)enrichGO()即可完成GO富集分析耻煤。
library(clusterProfiler)
#讀取基因列表文件中的基因名稱
genes <- read.delim('gene.txt', header = TRUE, stringsAsFactors = FALSE)[[1]]
#GO富集分析
#對(duì)于加載的注釋庫的使用买鸽，以上述為例，就直接在 OrgDb 中指定人（org.Hs.eg.db）或綿羊（sheep）
enrich.go <- enrichGO(gene = genes, #基因列表文件中的基因名稱
OrgDb = 'sheep', #指定物種的基因數(shù)據(jù)庫汽煮，示例物種是綿羊（sheep）
keyType = 'ENTREZID', #指定給定的基因名稱類型止后，例如這里以 entrze id 為例
ont = 'ALL', #可選 BP歉糜、MF、CC棵里，也可以指定 ALL 同時(shí)計(jì)算 3 者
pAdjustMethod = 'fdr', #指定 p 值校正方法
pvalueCutoff = 0.05, #指定 p 值閾值头谜，不顯著的值將不顯示在結(jié)果中
qvalueCutoff = 0.2, #指定 q 值閾值葵袭，不顯著的值將不顯示在結(jié)果中
readable = FALSE)
#例如上述指定 ALL 同時(shí)計(jì)算 BP鹉勒、MF、CC，這里將它們作個(gè)拆分后輸出
BP <- enrich.go[enrich.go$ONTOLOGY=='BP', ]
CC <- enrich.go[enrich.go$ONTOLOGY=='CC', ]
MF <- enrich.go[enrich.go$ONTOLOGY=='MF', ]
write.table(as.data.frame(BP), 'go.BP.txt', sep = '\t', row.names = FALSE, quote = FALSE)
write.table(as.data.frame(CC), 'go.CC.txt', sep = '\t', row.names = FALSE, quote = FALSE)
write.table(as.data.frame(MF), 'go.MF.txt', sep = '\t', row.names = FALSE, quote = FALSE)
對(duì)于各列內(nèi)容：
ONTOLOGY挤茄，GO分類BP（生物學(xué)過程）社证、CC（細(xì)胞組分）或MF（分子功能）疾渣；
ID和Description，富集到的GO id和描述项乒；
GeneRatio和BgRatio恋拍，分別為富集到該GO條目中的基因數(shù)目/給定基因的總數(shù)目僵娃，以及該條目中背景基因總數(shù)目/該物種所有已知的GO功能基因數(shù)目瘩扼；
pvalue、p.adjust和qvalue，p值仰猖、校正后p值和q值信息躏升；
geneID和Count卖局，富集到該GO條目中的基因名稱（分析中使用的entrze id染坯，故這里也顯示的entrze id）和數(shù)目湃鹊。
注：如期望顯示其它類型的基因id芯义，如通俗的symbol id等類型陶贼，除了更改為使用symbol id的基因名稱做分析外，還可以通過基因名稱轉(zhuǎn)換的方式對(duì)entrze id和symbol id作個(gè)匹配轉(zhuǎn)換培廓。
clusterProfiler的KEGG富集分析（有參向）
然后是KEGG富集分析当纱。
同樣地啥供，這里均對(duì)于KEGG數(shù)據(jù)庫中已經(jīng)收錄的物種而言的鳞骤，無參分析暫不涉及渤滞。
1唠粥、準(zhǔn)備輸入數(shù)據(jù)
相比上述GO富集只厘，clusterProfiler的KEGG富集分析方法特殊介评，它無需加載本地注釋庫垃你，自動(dòng)使用KEGG的在線數(shù)據(jù)庫進(jìn)行注釋，因此給定的基因名稱只能識(shí)別entrze id痴怨。把entrze id的基因名稱以一列的形式排開弓乙，放在一個(gè)文本文件中（例如命名“gene.txt”）乾颁。Excel中查看罚勾，就是如下示例這種樣式。
2、KEGG富集分析
讀取基因列表文件，并使用clusterProfiler的內(nèi)部函數(shù)enrichKEGG()即可完成KEGG富集分析跌榔。
library(clusterProfiler)
#讀取基因列表文件中的基因名稱示绊，注意這里只能用 entrze id
genes <- read.delim('gene.txt', header = TRUE, stringsAsFactors = FALSE)[[1]]
#每次打開R計(jì)算時(shí)摄杂，它會(huì)自動(dòng)連接kegg官網(wǎng)獲得最近的物種注釋信息映挂，因此數(shù)據(jù)庫一定都是最新的
kegg <- enrichKEGG(
gene = genes, #基因列表文件中的基因名稱
keyType = 'kegg', #kegg 富集
organism = 'oas', #例如亏拉，oas 代表綿羊灵再，其它物種更改這行即可
pAdjustMethod = 'fdr', #指定 p 值校正方法
pvalueCutoff = 0.05, #指定 p 值閾值，不顯著的值將不顯示在結(jié)果中
qvalueCutoff = 0.2, #指定 q 值閾值衷掷，不顯著的值將不顯示在結(jié)果中
)
#輸出結(jié)果
write.table(kegg, 'kegg.txt', sep = '\t', quote = FALSE, row.names = FALSE)
對(duì)于各列內(nèi)容：
ID和Description蚯根，富集到的KEGG id和描述；
GeneRatio和BgRatio，分別為富集到該KEGG條目中的基因數(shù)目/給定基因的總數(shù)目允趟，以及該條目中背景基因總數(shù)目/該物種所有已知的KEGG功能基因數(shù)目绽乔；
pvalue、p.adjust和qvalue，p值、校正后p值和q值信息逗余；
geneID和Count，富集到該KEGG條目中的基因名稱（分析中使用的entrze id，故這里也顯示的entrze id）和數(shù)目。
注：如期望顯示其它類型的基因id，如通俗的symbol id等類型，由于該分析中只能輸入entrze id菠齿，因此可以通過基因名稱轉(zhuǎn)換的方式對(duì)entrze id和symbol id作個(gè)匹配轉(zhuǎn)換疾棵。
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)系
