通過前面的講解，我們順利的了解了GEO數(shù)據(jù)庫以及如何下載其數(shù)據(jù)，得到我們想要的表達矩陣，也學(xué)會了兩個常用的套路分析得到的表達矩陣，就是GSEA分析和差異分析。
歷史目錄：
但是差異分析通過自定義的閾值挑選了有統(tǒng)計學(xué)顯著的基因列表后我們其實是需要對它們進行注釋才能了解其功能，最常見的就是GO/KEGG數(shù)據(jù)庫注釋咯，當(dāng)然也可以使用Reactome和Msigdb數(shù)據(jù)庫來進行注釋。而最常見的注釋方法就是超幾何分布檢驗了咯。
我2017年初我們生信技能樹舉辦的編程直播活動，我們專門講解了這個，看下面的第七題咯：
歷史題目：
ps：歷史題目都可以點擊閱讀原文可以查看其他同學(xué)提交的作業(yè)！
其實我在生信菜鳥團博客里面講的非常清楚了這個統(tǒng)計學(xué)原理：用超幾何分布檢驗做富集分析 http://www.bio-info-trainee.com/1225.html
直接上代碼： https://github.com/jmzeng1314/humanid/blob/master/R/hyperGtest_jimmy.R
維基百科的解釋是： 超幾何分布是統(tǒng)計學(xué)上一種離散概率分布。它描述了由有限個物件中抽出n個物件，成功抽出指定種類的物件的個數(shù)（不歸還）。
例如在有N個樣本，其中m個是不及格的。超幾何分布描述了在該N個樣本中抽出n個，其中k個是不及格的概率
若n=1，超幾何分布還原為伯努利分布。
若N接近∞，超幾何分布可視為二項分布。
作為離散概率分布的超幾何分布尤其指在抽樣試驗時抽出的樣品不再放回去的分布情況。在一個容器中一共有N個球，其中M個黑球,(N-M)個紅球，通過下面的超幾何分布公式可以計算出，從容器中抽出的n個球中（抽出的球不放回去）有k個黑球的概率是多少：
例如，容器中一共10個球，其中6個黑色，4個白色，一共抽5次（抽出的球不放回去），在這5個球中有3個黑球的概率是：
但是我們要算的不是恰好有3個黑球的概率，而是我們的基因富集問題，那么在R里面如何實現(xiàn)呢？
phyper(62-1, 1998, 5260-1998, 131, lower.tail=FALSE)
下面是我自己的理解：
超幾何分布很簡單，球分成黑白兩色，數(shù)量已知，那么你隨機抽有限個球，應(yīng)該抽多少白球的問題！
公式就是exp_count=n*M/N
然后你實際上抽了多少白球，就可以計算一個概率值！
換算成通路的富集概念就是，總共有多少基因（這個地方值得注意，主流認為只考慮那些在KEGG等數(shù)據(jù)庫注釋的背景基因），你的通路有多少基因，你的通路被抽中了多少基因（在差異基因里面屬于你的通路的基因），這樣的數(shù)據(jù)就足夠算出上面表格里面所有的數(shù)據(jù)啦！
原理講完了，就該講如何做了：https://github.com/jmzeng1314/humanid/blob/master/R/batch_enrichment.R
需要自己搜索什么是GO/KEGG/BIOCARTA/REACTOME等數(shù)據(jù)庫
http://www.cnblogs.com/emanlee/archive/2011/08/02/2125314.html
雖然懂了原理可以讓我們更方便的理解結(jié)果，但實際數(shù)據(jù)分析過程中我們通常不會自己寫代碼完成這個超幾何分布檢驗，因為有現(xiàn)成的，而且非常好用的輪子。
首先需要理解下面的 geneList和 gene這兩個數(shù)據(jù)集。然后，理解 GO/KEGG/REACTOME/MSIGDB 這4個數(shù)據(jù)庫結(jié)構(gòu)，及對應(yīng)的生物學(xué)一樣。接著，理解 超幾何分布建議，GSEA這兩個算法。最后把下面的代碼跑一遍即可。(因為Y叔的包更新太頻繁，我只能保證下面的代碼，今天是有效的，所以趕快試一下吧)
library(clusterProfiler)
### get the universal genes and sDEG
data(geneList, package='DOSE')
gene <- names(genelist)[abs(genelist)=""> 2]
gene.df <- bitr(gene,="" fromtype='ENTREZID'>
toType = c('ENSEMBL', 'SYMBOL'),
OrgDb = org.Hs.eg.db)
head(gene.df)
### Using MSigDB gene set collections
gmtfile <- system.file('extdata',="" 'c5.cc.v5.0.entrez.gmt',="" package='clusterProfiler'>
c5 <->
egmt2 <- gsea(genelist,="" term2gene="c5," verbose="">
head(egmt2)[,1:6]
gseaplot(egmt2, geneSetID = 'EXTRACELLULAR_REGION')
## GO analysis
ego <- enrichgo(gene="" ="" ="" ="" ="" ="">
universe = names(geneList),
OrgDb = org.Hs.eg.db,
ont = 'CC',
pAdjustMethod = 'BH',
pvalueCutoff = 0.01,
qvalueCutoff = 0.05,
readable = TRUE)
head(ego)[,1:6]
ego2 <- gsego(genelist="" ="" ="">
OrgDb = org.Hs.eg.db,
ont = 'CC',
nPerm = 1000,
minGSSize = 100,
maxGSSize = 500,
pvalueCutoff = 0.05,
verbose = FALSE)
head(ego2)[,1:6]
## KEGG pathway analysis
kk <- enrichkegg(gene="" ="" ="" ="" ="">
organism = 'hsa',
pvalueCutoff = 0.05)
head(kk)[,1:6]
kk2 <- gsekegg(genelist="" ="" ="">
organism = 'hsa',
nPerm = 1000,
minGSSize = 120,
pvalueCutoff = 0.05,
verbose = FALSE)
head(kk2)[,1:6]
gseaplot(kk2, geneSetID = 'hsa04145')
### Reactome pathway analysis
library(ReactomePA)
pp <- enrichpathway(gene="" ="" ="" ="" ="">
organism = 'human',
pvalueCutoff = 0.05)
head(pp)[,1:6]
pp2 <- gsepathway(genelist="" ="" ="">
organism = 'human',
nPerm = 1000,
minGSSize = 120,
pvalueCutoff = 0.05,
verbose = FALSE)
head(pp2)[,1:6]
### Disease analysis
library(DOSE)
dd <- enrichdo(gene="" ="" ="" ="" ="">
head(dd)[,1:6]
dd2 <- gsedo(genelist="" ="" ="">
nPerm = 1000,
minGSSize = 120,
pvalueCutoff = 0.05,
verbose = FALSE)
head(dd2)[,1:6]
下游的GO/KEGG注釋一般是得到如下表格：
也可以有一些簡單的可視化展現(xiàn)：
