我們統(tǒng)一選擇p<0.05而且abs(logFC)大于1的基因?yàn)轱@著差異表達(dá)基因集前联，對這個基因集用R包做KEGG/GO超幾何分布檢驗(yàn)分析支鸡。
然后把表達(dá)矩陣和分組信息分別作出cls和gct文件悦陋，導(dǎo)入到GSEA軟件分析慕爬。
基本任務(wù)是完成這個分析
其實(shí)這一步特別的簡單畅铭，就是篩選淡溯，然后用專門的R包分析就好了。但是富集分析貌似簡單拒逮，但其實(shí)充滿了變數(shù)罐氨。
PS: 下面的內(nèi)容我直接從我之前的文章里摘錄過來了。
為何要基因富集分析
在基因差異表達(dá)分析之后滩援，你得到了好多p值特別姓ひ（也就是顯著性很高）的基因，那么下一步你想做什么？
- 選擇一些基因用于驗(yàn)證租悄？
- 對其中基因進(jìn)行后續(xù)研究谨究？
- 在結(jié)果中把這些基因都放在后面？
- 嘗試著把所有基因相關(guān)的文獻(xiàn)都都讀讀看（勸你放棄這個念頭）泣棋？
- 歡迎補(bǔ)充
這些想法都是非常順理成章的胶哲，但是不要著急。
首先潭辈，差異表達(dá)找到的基因往往很多鸯屿，你簡單的粗暴去找每一個基因的詳細(xì)資料，顯然不太現(xiàn)實(shí)把敢；
其次寄摆，如果我們單純覺得某一個基因和你研究的課題相關(guān)，或者說你其實(shí)已經(jīng)找到了一個有可能的基因（或者你只是希望用一些高大上的實(shí)驗(yàn)驗(yàn)證一下）那么這個行為是不是有太多主觀性修赞，存在一些偏見婶恼。
當(dāng)然，你覺得基因就是你要找的柏副，可是萬一它只是碰巧來打醬油的呢勾邦，這不就是很尷尬了。
所以為了讓審稿人相信你的結(jié)果搓扯，你就需要做一個基因富集分析哦检痰。
什么是基因富集分析
基因富集分析(gene set enrichment analysis)是在一組基因或蛋白中找到一類過表達(dá)的基因或蛋白包归。一般是高通量實(shí)驗(yàn)锨推，如基因芯片，RNA-Seq公壤，蛋白質(zhì)組學(xué)（質(zhì)譜結(jié)果）的后續(xù)步驟换可。
基因富集分析需要我們提供某一類功能基因的集合用于背景，常用的注釋數(shù)據(jù)庫如：
- The Gene Ontology Consortium: 描述基因的層級關(guān)系
- Kyoto Encyclopedia of Genes and Genomes: 提供了pathway的數(shù)據(jù)庫厦幅。
分析方法
在文獻(xiàn)Ten Years of Pathway Analysis: Current Approaches and Outstanding Challenges（推薦大家看一遍）作者將研究方法歸為三種沾鳄，其中第三種方法想的很好就是難度很大。：
作者還貼心的把每一種方法有哪些工具都總結(jié)出來了：
Over-Repressentation Analysis(ORA)
ORA是目前商業(yè)化最多的方法确憨。為了說明他的基本思想译荞，我要舉一個喜聞樂見的例子：讀書無用論。
這是我百度找到搜狐財(cái)經(jīng)一篇文章《大數(shù)據(jù)告訴你真正的有錢人是什么樣》的有錢人的學(xué)歷分布情況休弃，高學(xué)歷人群（本科以上吞歼，因?yàn)楸究粕嗔耍┧嫉谋壤?.4%，其他都是一般學(xué)歷占90.6%塔猾。這時候篙骡，有些公眾號就可以開始不帶腦子的說了，讀書沒什么用呀，有錢人中都是一般學(xué)歷的呀糯俗，以后讀書讀到大學(xué)就行了尿褪，甚至也可以不上本科呀（34.2%本科和本科以下）。
你每年回家總能回去看到有人炫耀說，雖然我有錢崭歧，可是讀書太少了诞帐，都不能和你們讀書人比的。你總感覺哪里不對勁天揖，但是卻又不太方便說出來。
實(shí)際上跪帝，這就是因?yàn)闆]有考慮到背景今膊。因?yàn)楦邔W(xué)歷本身人數(shù)就不多，當(dāng)然在有錢人里面的人數(shù)也就相應(yīng)不多了伞剑。我們要證明有錢人更多是富集高學(xué)歷這一部分斑唬。
| 類別 | 有錢人 | 普通人 |
|---|---|---|
| 高學(xué)歷 | 10 | 50 |
| 一般學(xué)歷 | 90 | 850 |
| 100 | 900 |
H0: 是否有錢和學(xué)歷高無關(guān)
Ha: 學(xué)歷高還是有點(diǎn)用的
然后做一個Fisher精確檢驗(yàn)，看看p值黎泣。
richer.pop <- matrix(data = c(10,90,50,850),nrow=2)
fisher.test(richer.pop, alternative = "greater")
Fisher's Exact Test for Count Data
data: richer.pop
p-value = 0.03857
alternative hypothesis: true odds ratio is greater than 1
95 percent confidence interval:
1.052584 Inf
sample estimates:
odds ratio
2.109244
p值小于0.05恕刘，看來我讀個博士讓我以后有錢概率變大了。
現(xiàn)在將我們上面的有錢人改成我們找到的基因抒倚，整體改成所有基因褐着。高學(xué)歷表示屬于目標(biāo)注釋基因集，一般學(xué)歷就是非注釋基因組.我們就是要判斷我們找到的基因更多是在目標(biāo)注釋集中托呕。所以你需要列出下表含蓉，然后再做一個fisher.test()。
| 類別 | 感興趣的基因 | 其他基因 |
|---|---|---|
| in anno group | 10 | 50 |
| not in anno group | 290 | 19950 |
| 300 | 20000 |
上述的基本思想就是統(tǒng)計(jì)學(xué)的白球黑球?qū)嶒?yàn)：
在一個黑箱里项郊，有確定數(shù)量的黑白兩種球馅扣，你隨機(jī)抽取（不放回）M個球中着降，其中兩種球的比例分別是多少差油？
除了用Fisher精確檢驗(yàn)，還有其他統(tǒng)計(jì)方法：
- Hypergeometric (fisher精確檢驗(yàn)用的就是超幾何檢驗(yàn))http://www.bio-info-trainee.com/1225.html
- Binomial: 二項(xiàng)分布要求是有放回任洞，無放回要求整體足夠大大到可以近似蓄喇。
- Chi-squared
chisq.test(counts) - Z
- Kolmogorov-Smirnov
- Permutation http://www.bio-info-trainee.com/1237.html
ORA的方法就是如此的簡單，但是有一個問題交掏，就是你如何確定哪些基因是差異表達(dá)的妆偏，你還是需要設(shè)置一個人為的cutoff, 主觀能動性成分有點(diǎn)大。
Functional Class Scoring(FCS)
FCS認(rèn)為耀销，“雖然個體基因表達(dá)改變之后會更多在通路中體現(xiàn)楼眷，但是一些功能相關(guān)基因中較弱但協(xié)調(diào)的變化也有明顯的影響铲汪。”
The hypothesis of functional class scoring (FCS) is that although large changes in individual genes can have significant effects on pathways, weaker but coordinated changes in sets of functionally related genes (i.e., pathways) can also have significant effects
FCS分析方法稍微復(fù)雜了一點(diǎn)罐柳，他要求的輸入是一個排序的基因列表和一個基因集合掌腰。MIT， Broad Institute 2007年文獻(xiàn)就提供了這一方法的軟件"GSEA"
有如下特點(diǎn)：
- 計(jì)算所有輸入基因集合的分?jǐn)?shù)张吉，而不是單個基因
- 不需要設(shè)置cutoff
- 找到一組相關(guān)的基因
- 提供了更加穩(wěn)健的統(tǒng)計(jì)框架
GSEA是一款圖形化的軟件齿梁，根據(jù)他們提供的教程，然后點(diǎn)呀點(diǎn)肮蛹，就會得到如下結(jié)果勺择。下圖就是需要好好理解的部分。
中間從藍(lán)色到紅色的過渡“帶”表示基因從上調(diào)到下調(diào)排列（排序可以按照fold change,也可以是p-value)伦忠。黑色像條形碼的豎線表示該位置的基因?qū)儆谀硞€指定通路省核。綠色有波動的曲線表示富集分?jǐn)?shù)，從0開始計(jì)算昆码，屬于基因通路增加气忠，不屬于則減少。最后看下黑色的條形碼是不是富集在一端赋咽。
那如何做統(tǒng)計(jì)檢驗(yàn)?zāi)兀?/p>
The final step in FCS is assessing the statistical significance of the pathway-level statistic. When computing statistical significance, the null hypothesis tested by current pathway analysis approaches can be broadly divided into two categories: i) competitive null hypothesis and ii) self-contained null hypothesis [3], [18], [22], [31]. A self-contained null hypothesis permutes class labels (i.e., phenotypes) for each sample and compares the set of genes in a given pathway with itself, while ignoring the genes that are not in the pathway. On the other hand, a competitive null hypothesis permutes gene labels for each pathway, and compares the set of genes in the pathway with a set of genes that are not in the pathway旧噪。
我們要檢驗(yàn)的目標(biāo)是基因富集在一端是因?yàn)橛谀繕?biāo)通路相關(guān)的基因都在一端富集。那么空假設(shè)就是脓匿，你把找到的基因隨便擺放也能看到富集現(xiàn)象淘钟。用比較專業(yè)的話說就是先生成一個零假設(shè)的數(shù)據(jù)分布，然后觀察實(shí)際數(shù)據(jù)在這個零假設(shè)分布下陪毡，是不是在尾端 米母。
一些問題
統(tǒng)計(jì)檢驗(yàn)的能力是有限的，所以還有很多問題存在解決缤骨。
- 我們希望找到生物學(xué)顯著的基因爱咬，但是生物學(xué)顯著和統(tǒng)計(jì)顯著兩者并不是完全相關(guān)
- 無論是ORA還是FCS都對背景（也就是這個物種一共有多少基因）有要求，但是隨著我們的研究深入绊起，基因數(shù)量會改變。有些軟件會直接設(shè)置一個很大的背景數(shù)燎斩，從而讓p值很顯著虱歪，然后我們就開心地用他們的結(jié)果。
- 有些基因沒有注釋栅表，也就是注釋缺失笋鄙，處理方法就是扔（歡迎拍磚）。
- 有一些注釋項(xiàng)是其他項(xiàng)的子集怪瓶。
富集分析
我們這里用來做富集分析的工具是Y叔的clusterProfiler萧落。基本上用這個軟件就對了，簡單還用找岖，而且結(jié)果正確陨倡。它支持上面提到的三類方法，所以也沒有必要專門導(dǎo)出數(shù)據(jù)用GSEA軟件了：
- Over-Representation Analysis
- Gene Set Enrichment Analysis
- Biological theme comparison
前提準(zhǔn)備：
- 數(shù)據(jù)篩選许布，根據(jù)padj < 0.05 且Log2FoldChange的絕對值大于1的標(biāo)準(zhǔn)兴革。
deseq2.sig <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
- 安裝包下載注釋數(shù)據(jù)（rog.HS.eg.db)
source("https://bioconductor.org/biocLite.R")
biocLite("clusterProfiler")
biocLite（"AnnotationHub")
library(AnnotationHub)
ah <- AnnotationHub()
org.hs <- ah[['AH53766']]
注釋數(shù)據(jù)庫根據(jù)分析的物種決定，如何查找數(shù)據(jù)庫蜜唾，可以看我寫的一篇文章用Bioconductor對基因組注釋# 用Bioconductor對基因組注釋
GO富集和GESA
GO分析參考Y叔寫的GO analysis using clusterProfiler文檔：
主要的函數(shù)是
enrichGO(gene, OrgDb, keytype = "ENTREZID", ont = "MF",
pvalueCutoff = 0.05, pAdjustMethod = "BH", universe, qvalueCutoff = 0.2,
minGSSize = 10, maxGSSize = 500, readable = FALSE, pool = FALSE)
主要關(guān)注如下參數(shù)
gene: 差異表達(dá)的基因. 不推薦使用"A1BG"這中類型的命名杂曲，請用setReadable轉(zhuǎn)換
OrgDb： 物種注釋數(shù)據(jù)庫，一般是org開頭
keytpe: gene的命名格式
ont： 是BP（Biological Process）, CC（Cellular Component）, MF（Molecular Function）袁余。一個基因的功能可以從生物學(xué)過程擎勘，所屬細(xì)胞部分，和分子功能三個角度定義
只要提供對應(yīng)的數(shù)據(jù)即可
ego <- enrichGO(
gene = row.names(deseq2.sig),
OrgDb = org.hs,
keytype = "ENSEMBL",
ont = "MF"
)
可視化分為冒泡圖和網(wǎng)絡(luò)圖等颖榜，畫起來也就幾行代碼的事情
dotplot(ego,font.size=5)
enrichMap(ego, vertex.label.cex=1.2, layout=igraph::layout.kamada.kawai)
plotGOgraph(ego)
<p align="center">氣泡圖</p>
<p align="center">網(wǎng)絡(luò)圖</p>
<p align="center">GO圖</p>
GSEA分析可使用broadinstitute出品的GSEA可視化軟件包.有興趣看下Y叔的文章Comparison of clusterProfiler and GSEA-P
因此我們用clusterProfiler的gseGO或GSEA函數(shù)分析货抄，后者可以自定義輸入數(shù)據(jù)
gseGO(geneList, ont = "BP", OrgDb, keyType = "ENTREZID", exponent = 1,
nPerm = 1000, minGSSize = 10, maxGSSize = 500, pvalueCutoff = 0.05,
pAdjustMethod = "BH", verbose = TRUE, seed = FALSE, by = "fgsea")
geneList： 排序數(shù)據(jù)， 可以根據(jù)log2foldchange, 也可以是pvalues
nPerm： 重抽取次數(shù)
minGSSize： 每個基因集的最小數(shù)目
maxGSSize： 用于測試的基因注釋最大數(shù)目
genelist <- sig.deseq2$log2FoldChange
names(genelist) <- rownames(sig.deseq2)
genelist <- sort(genelist, decreasing = TRUE)
gsemf <- gseGO(genelist,
OrgDb = org.hs,
keyType = "ENSEMBL",
ont="MF"
)
head(gsemf)
畫一個GSEA標(biāo)志性圖
gseaplot(gsemf, geneSetID="GO:0004871")
<p align="center">GO:0004871</p>
KEGG富集分析
KEGG富集分析那家強(qiáng)朱转，必須是Y叔的clusterProfiler蟹地。因?yàn)樗芘廊∽钚碌腒EGG在線版數(shù)據(jù)庫，而不是用不再更新的KEGG.db藤为。
本部分內(nèi)容參照KEGG enrichment analysis with latest online data using clusterProfiler
函數(shù)是enrichKEGG
enrichKEGG(gene, organism = "hsa", keyType = "kegg", pvalueCutoff = 0.05,
pAdjustMethod = "BH", universe, minGSSize = 10, maxGSSize = 500,
qvalueCutoff = 0.2, use_internal_data = FALSE)
- gene： 基因名怪与，要和keyType對應(yīng)
- organism: 需要參考 http://www.genome.jp/kegg/catalog/org_list.html， 人類是hsa
- keyType: 基因的命名方式缅疟， "kegg", 'ncbi-geneid', 'ncib-proteinid' and 'uniprot'選擇其一
library(clusterProfiler)
gene_list <- mapIds(org.hs, keys = row.names(deseq2.sig),
column = "ENTREZID", keytype = "ENSEMBL" )
kk <- enrichKEGG(gene_list, organism="hsa",
keyType = "ncbi-geneid",
pvalueCutoff=0.05, pAdjustMethod="BH",
qvalueCutoff=0.1)
head(summary(kk))
# 氣泡圖等圖限于篇幅不畫了分别，流量有限
dotplot(kk)
注意： 一定要提供符合keyType要求的基因名
總結(jié)
你以為結(jié)束了？這只是剛開始存淫！我還得修修補(bǔ)補(bǔ)耘斩，說不定能出書呢！
參考文獻(xiàn)
[1] Guangchuang Yu., Li-Gen Wang, Yanyan Han, Qing-Yu He. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS: A Journal of Integrative Biology. 2012, 16(5):284-287.
[2] Ten Years of Pathway Analysis: Current Approaches and Outstanding Challenges
