1欠母、基礎(chǔ)知識(shí)
（1）基本概念
富集分析(enrichment analysis)簡單來說就是將成百上千個(gè)基因皆辽、蛋白或者其他分子分到不同的類中莹菱，以減少分析的復(fù)雜度。比如之前差異分析得到幾百個(gè)顯著差異基因蛤育，如果一個(gè)一個(gè)單獨(dú)研究未免太復(fù)雜，若按照一定的準(zhǔn)則將差異基因歸類即可較為快速葫松，方便的了解某一類基因的變化情況瓦糕。
（2）分類標(biāo)準(zhǔn)
分類的標(biāo)準(zhǔn)即人們根據(jù)目前研究建立的基因注釋庫，目前常用的有兩個(gè)：GO與KEGG腋么；
- GO簡單來說GO term共有三種類型 ①細(xì)胞組成（cellular component咕娄，CC）；②生物過程（biological process党晋，BP）谭胚；③分子功能（Molecular Function，MF）未玻。每個(gè)go term都由相應(yīng)的GO annotation來說明該term的詳細(xì)信息灾而，例如人的go注釋文件為org.Hs.eg.db（該類型文件均為.db結(jié)尾）。通過GO富集分析可以了解差異基因富集在哪些生物學(xué)功能扳剿、途徑或細(xì)胞定位旁趟。
- KEGG，Kyoto Encyclopedia of Genes and Genomes庇绽，京都基因和基因組百科全書锡搜，是系統(tǒng)分析基因功能，聯(lián)系基因組信息和功能信息的知識(shí)庫,其中包含有大量的通路（PATHWAY）圖瞧掺。人的KEGG注釋文件為“hsa”耕餐。KEGG分析的最終結(jié)果就是把判斷某些基因是都富集到某一通路上。
（3）過表達(dá)分析 over-representation analysis
舉一例：首先獲得一組感興趣的基因(一般是差異表達(dá)基因辟狈，前景基因)肠缔，假設(shè)有10個(gè)；其中有4個(gè)都?xì)w類到某一GO term中或者落在某一通路中（pathway）哼转；而在整個(gè)基因組中（假設(shè)為100個(gè)明未，背景基因）有30個(gè)都對(duì)應(yīng)該GO term中或者落在該通路中（pathway）中∫悸基于此來研究4/10與30/100間是否有統(tǒng)計(jì)學(xué)差異趟妥，即觀察的計(jì)數(shù)值是否顯著高于隨機(jī)，即待測功能集在基因列表中是否顯著富集佣蓉。
（4）分析平臺(tái)
目前有蠻多不錯(cuò)的網(wǎng)站在線富集分析軟件披摄，當(dāng)然也可以通過R語言的R包實(shí)現(xiàn)亲雪。這里以眾多人推薦的clusterProfiler包為例進(jìn)行學(xué)習(xí)
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("clusterProfiler")
library(clusterProfiler)
差異基因數(shù)據(jù)使用前期基于airway包分析得到的結(jié)果(篩選條件為
padj<0.1&abs(log2FoldChange)>2)
mydata=read.table("results.csv",header=TRUE,
sep=",",stringsAsFactors=FALSE)
genelist=mydata$X
genelist #共有316個(gè)差異基因
注意此時(shí)的基因名為
ENSEMBL格式，特征是以ENSG00000字段開頭的行疏；其它常見的格式還有ENTREZID匆光，為純數(shù)字序列；SYMBOL為字母為主的字符串酿联。
2终息、go分析
BiocManager::install("org.Hs.eg.db")
library("org.Hs.eg.db")
# 安裝加載人類go注釋包
go.all <- enrichGO(genelist, OrgDb = org.Hs.eg.db, ont='ALL',pAdjustMethod = 'BH',pvalueCutoff = 0.05,
qvalueCutoff = 0.2,keyType = 'ENSEMBL')
# 需要稍等一會(huì)
-
enrichGO函數(shù)可以支持多種基因名格式，使用keyType =指定下即可贞让。
dim(go.all)
head(go.all)
go.all.df=as.data.frame(go.all)
由返回結(jié)果可以看出共歸類到97個(gè) go term中周崭，其中BP類較多；關(guān)于表格中的部分變量意義--
-
Description：Gene Ontology功能的描述信息; -
GeneRatio：差異基因中與該Term相關(guān)的基因數(shù)與整個(gè)差異基因總數(shù)的比值喳张； -
BgRation：所有背景基因中與該Term相關(guān)的基因數(shù)與所有基因的比值续镇； - 三個(gè)value值：一般情況下， pvalue < 0.05 該功能為富集項(xiàng)销部；p.adjust 矯正后的pvalue摸航；qvalue為對(duì)p值進(jìn)行統(tǒng)計(jì)學(xué)檢驗(yàn)的q值；
-
Count：差異基因中與該Term相關(guān)的基因數(shù)舅桩。
go.all.df
genelist里的差異基因明明有316個(gè)酱虎，表格中顯示好像只有221個(gè)？擂涛？奇怪
基于上述分析读串，將結(jié)果可視化
（1）柱形圖
barplot(go.all,showCategory=20,drop=T)
- 篩選20個(gè)p.adjust值最小的GO term；
-
統(tǒng)計(jì)量為Count值撒妈；顏色程度根據(jù)p.adjust恢暖；縱軸標(biāo)簽為Description
barplot
（2）點(diǎn)圖
dotplot(go.aal,showCategory=20)
- 篩選20個(gè)分類到某一term中基因數(shù)最多的GO term
-
注意下GeneRatio是與Count值成正比的，可查看之前的定義
dotplot
（3）有向無環(huán)圖DAG
由于這里的數(shù)據(jù)只能是富集一個(gè)GO通路（BP狰右、CC或MF）的數(shù)據(jù)杰捂，因此重新針對(duì)某一類go，再分析一下棋蚌。
go.BP <- enrichGO(genelist,
OrgDb = org.Hs.eg.db,
ont='BP',
pAdjustMethod = 'BH',
pvalueCutoff = 0.05,
qvalueCutoff = 0.2,
keyType = 'ENSEMBL')
plotGOgraph(go.BP)
# 之前失敗琼娘，提示說一下兩個(gè)相關(guān)包未找到
# BiocManager::install("topGO")
# BiocManager::install("Rgraphviz")
-
得到的結(jié)果圖比較大，建議保存成pdf文件查看
plotGOgraph(go.BP)
#覺得默認(rèn)文字有點(diǎn)小附鸽，調(diào)大一下~
opar <- par(no.readonly=TRUE)
par.axis(cex=4)
plotGOgraph(go.BP)
par(opar)
如上圖DAG圖結(jié)果--
- 顏色越深，代表該GO term越顯著（p值越新魅场）坷备，函數(shù)默認(rèn)將最顯著的10個(gè)term設(shè)置成方形；
- 圖形內(nèi)標(biāo)注信息分別是GOterm號(hào)情臭、Description省撑、P.adjust以及差異基因注釋到該term數(shù)與背景基因注釋到該term數(shù)的比值赌蔑；
- 超接近根結(jié)點(diǎn)的GO term的概括性越強(qiáng)，越往下竟秫，分支的GO term表示的結(jié)果更細(xì)娃惯。
3、kegg分析
由于enrichKEGG()需要輸入的基因名格式為ENTREZID肥败，所以需要轉(zhuǎn)換一下趾浅，這里使用clusterProfiler包的bitr()函數(shù)
gene=bitr(genelist,fromType="ENSEMBL", toType="ENTREZID", OrgDb="org.Hs.eg.db")
gene=gene$ENTREZID
基因數(shù)有316變成了267，沒匹配到馒稍？
kegg <- enrichKEGG(gene,
organism = 'hsa',
keyType = 'kegg',
pvalueCutoff = 0.05,
pAdjustMethod = 'BH',
minGSSize = 10,
maxGSSize = 500,
qvalueCutoff = 0.2
)
head(kegg)
dim(kegg)
kegg.df=as.data.frame(kegg)
-
結(jié)果好像只有富集到兩條通路上皿哨，由于第一次做，不知道對(duì)于我的數(shù)據(jù)來說結(jié)果是否正常纽谒。
kegg.df
列的意義可參考go分析得到的表格
結(jié)果可視化
（1）柱狀圖证膨、點(diǎn)圖
# 畫法同前
barplot(kegg)
dotplot(kegg)
（2）通路圖--針對(duì)每一個(gè)富集到的通路圖畫的
browseKEGG(kegg,'hsa04213')
-
紅色邊框的為上調(diào)的，綠色邊框的為下調(diào)的鼓黔。
部分截圖
以上就是基于ORP方法央勒，利用兩種注釋庫（go、kegg）進(jìn)行的分析澳化。其實(shí)過表達(dá)分析有不少缺點(diǎn)的崔步，比如
（1）僅使用了基因數(shù)目信息，而沒有利用基因表達(dá)水平或表達(dá)差異值肆捕；而基因篩選條件基于人為選定差異水平（比如log2FoldChange）刷晋；
（2）忽略差異不顯著，但比較關(guān)鍵的基因慎陵；
（3）將基因同等對(duì)待眼虱，ORA法假設(shè)每個(gè)基因都是獨(dú)立的，忽視了基因在通路內(nèi)部生物學(xué)意義的不同(如調(diào)控和被調(diào)控基因的不同)及基因間復(fù)雜的相互作用席纽；
（4）ORA假設(shè)通路與通路間是獨(dú)立的捏悬，但這個(gè)前提假設(shè)是錯(cuò)誤的。由于上述的不足润梯，GSEA方法也常為人們的選擇过牙，將會(huì)在下一次中繼續(xù)總結(jié)。
參考文章
1纺铭、功能富集分析概述
2寇钉、GO分析學(xué)習(xí)筆記（推薦!）
