回顧
首先對(duì)基因表達(dá)分析（上）做一個(gè)簡單的回顧
- 研究基因表達(dá)的有如下工具：RNA-Seq峡捡，microarray擦酌， qRT-PCR等（歡迎補(bǔ)充）
- RNA-Seq翁脆，microarray一般用在探索性階段筹裕，qRT-PCR用于驗(yàn)證
- RNA-Seq和microarray由于他們的實(shí)驗(yàn)方式不同抖僵，導(dǎo)致尋找差異表達(dá)基因的統(tǒng)計(jì)學(xué)方法也不同。其中microarray使用寡核苷酸作為探針進(jìn)行雜交宅广，基因表達(dá)量與亮度正相關(guān)葫掉，而亮度是一個(gè)連續(xù)型變量，因此大多認(rèn)為結(jié)果是服從正態(tài)分布跟狱。而RNA-Seq的測(cè)序結(jié)果是一條條read，是一種離散抽樣過程户魏，因此認(rèn)為是服從泊松分布驶臊。
- ANOVA和簡單線性模型都是廣義線性模型的特殊情況。ANOVA是研究名義型解釋變量和連續(xù)型解釋變量的關(guān)系叼丑，簡單線性模式是研究連續(xù)型解釋變量和連續(xù)型解釋變量的關(guān)系关翎。而廣義線性模式?jīng)]特殊要求。
- 在3,4的背景下鸠信，microarray一般用t檢驗(yàn)（兩個(gè)條件）纵寝，ANOVA分析（多個(gè)條件），最常用limma（線性模型）進(jìn)行檢驗(yàn)星立。RNA-Seq有許多基于count的R包爽茴，如DESeq，DESeq2,（基于負(fù)二向分布廣義線性模型）
- 以上要求你每個(gè)條件都要有3個(gè)重復(fù)（目前投稿要求）绰垂，你要是老板窮室奏，一個(gè)重復(fù)都不給，那你去Google解決方案吧劲装。
- 用R作差異表達(dá)分析大致分為以下幾步：1）根據(jù)軟件包要求導(dǎo)入數(shù)據(jù)胧沫；2）數(shù)據(jù)預(yù)處理昌简，把那些只有0或1計(jì)數(shù)結(jié)果的基因去掉，提高效率绒怨。這一步還可以進(jìn)行探索性數(shù)據(jù)分析纯赎；3）跑程序，得到結(jié)果南蹂；4）對(duì)結(jié)果進(jìn)行可視化犬金，看看基因聚類等結(jié)果，這一步不是必須的碎紊，但卻是展示數(shù)據(jù)最好的手段了佑附。
如果這些內(nèi)容都還有印象，那么我們繼續(xù)上次沒有寫完的內(nèi)容的其中一個(gè)部分--基因富集分析
為何要基因富集分析
在基因差異表達(dá)分析之后仗考，你得到了好多p值特別幸敉（也就是顯著性很高）的基因，那么下一步你想做什么秃嗜？
- 選擇一些基因用于驗(yàn)證权均？
- 對(duì)其中基因進(jìn)行后續(xù)研究？
- 在結(jié)果中把這些基因都放在后面锅锨？
- 嘗試著把所有基因相關(guān)的文獻(xiàn)都都讀讀看（勸你放棄這個(gè)念頭）叽赊？
- 歡迎補(bǔ)充
這些想法都是非常順理成章的，但是不要著急必搞。
首先必指，差異表達(dá)找到的基因往往很多，你簡單的粗暴去找每一個(gè)基因的詳細(xì)資料恕洲，顯然不太現(xiàn)實(shí)塔橡；
其次，如果我們單純覺得某一個(gè)基因和你研究的課題相關(guān)霜第，或者說你其實(shí)已經(jīng)找到了一個(gè)有可能的基因（或者你只是希望用一些高大上的實(shí)驗(yàn)驗(yàn)證一下）那么這個(gè)行為是不是有太多主觀性葛家，存在一些偏見。
當(dāng)然泌类，你覺得基因就是你要找的癞谒，可是萬一它只是碰巧來打醬油的呢，這不就是很尷尬了刃榨。
所以為了讓審稿人相信你的結(jié)果弹砚，你就需要做一個(gè)基因富集分析哦。
什么是基因富集分析
基因富集分析(gene set enrichment analysis)是在一組基因或蛋白中找到一類過表達(dá)的基因或蛋白喇澡。一般是高通量實(shí)驗(yàn)迅栅，如基因芯片，RNA-Seq晴玖，蛋白質(zhì)組學(xué)（質(zhì)譜結(jié)果）的后續(xù)步驟读存。
基因富集分析需要我們提供某一類功能基因的集合用于背景为流，常用的注釋數(shù)據(jù)庫如：
- The Gene Ontology Consortium: 描述基因的層級(jí)關(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è)化最多的方法敬察。為了說明他的基本思想，我要舉一個(gè)喜聞樂見的例子：讀書無用論尔当。
這是我百度找到搜狐財(cái)經(jīng)一篇文章《大數(shù)據(jù)告訴你真正的有錢人是什么樣》的有錢人的學(xué)歷分布情況莲祸，高學(xué)歷人群（本科以上，因?yàn)楸究粕嗔耍┧嫉谋壤?.4%椭迎，其他都是一般學(xué)歷占90.6%揩局。這時(shí)候断傲，有些公眾號(hào)就可以開始不帶腦子的說了贴妻，讀書沒什么用呀簇抵，有錢人中都是一般學(xué)歷的呀，以后讀書讀到大學(xué)就行了简软，甚至也可以不上本科呀（34.2%本科和本科以下）蛮拔。
你每年回家總能回去看到有人炫耀說，雖然我有錢痹升，可是讀書太少了建炫，都不能和你們讀書人比的。你總感覺哪里不對(duì)勁疼蛾，但是卻又不太方便說出來肛跌。
實(shí)際上，這就是因?yàn)闆]有考慮到背景察郁。因?yàn)楦邔W(xué)歷本身人數(shù)就不多惋砂，當(dāng)然在有錢人里面的人數(shù)也就相應(yīng)不多了。我們要證明有錢人更多是富集高學(xué)歷這一部分绳锅。
| 類別 | 有錢人 | 整體 |
|---|---|---|
| 高學(xué)歷 | 10 | 50 |
| 一般學(xué)歷 | 90 | 950 |
| 100 | 1000 |
H0: 是否有錢和學(xué)歷高無關(guān)
Ha: 學(xué)歷高還是有點(diǎn)用的
然后做一個(gè)Fisher精確檢驗(yàn)，看看p值酝掩。
richer.pop <- matrix(data = c(10,90,50,950),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鳞芙，看來我讀個(gè)博士讓我以后有錢概率變大了。
現(xiàn)在將我們上面的有錢人改成我們找到的基因期虾，整體改成所有基因原朝。高學(xué)歷表示屬于目標(biāo)注釋基因集，一般學(xué)歷就是非注釋基因組.我們就是要判斷我們找到的基因更多是在目標(biāo)注釋集中镶苞。所以你需要列出下表喳坠，然后再做一個(gè)fisher.test()。
| 類別 | gene list | Genome |
|---|---|---|
| in anno group | 10 | 50 |
| not in anno group | 290 | 19950 |
| 300 | 20000 |
上述的基本思想就是統(tǒng)計(jì)學(xué)的白球黑球?qū)嶒?yàn)：
在一個(gè)黑箱里茂蚓，有確定數(shù)量的黑白兩種球壕鹉，你隨機(jī)抽忍昊稀（不放回）M個(gè)球中，其中兩種球的比例分別是多少晾浴？
除了用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的方法就是如此的簡單抖棘，但是有一個(gè)問題，就是你如何確定哪些基因是差異表達(dá)的狸涌，你還是需要設(shè)置一個(gè)人為的cutoff, 主觀能動(dòng)性成分有點(diǎn)大切省。
Functional Class Scoring(FCS)
FCS認(rèn)為，“雖然個(gè)體基因表達(dá)改變之后會(huì)更多在通路中體現(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)惶楼，他要求的輸入是一個(gè)排序的基因列表和一個(gè)基因集合右蹦。MIT， Broad Institute 2007年文獻(xiàn)就提供了這一方法的軟件"GSEA"
有如下特點(diǎn)：
- 計(jì)算所有輸入基因集合的分?jǐn)?shù)歼捐，而不是單個(gè)基因
- 不需要設(shè)置cutoff
- 找到一組相關(guān)的基因
- 提供了更加穩(wěn)健的統(tǒng)計(jì)框架
GSEA是一款圖形化的軟件何陆，根據(jù)他們提供的教程，然后點(diǎn)呀點(diǎn)豹储，就會(huì)得到如下結(jié)果贷盲。下圖就是需要好好理解的部分。
中間從藍(lán)色到紅色的過渡“帶”表示基因從上調(diào)到下調(diào)排列（排序可以按照fold change,也可以是p-value)剥扣。黑色像條形碼的豎線表示該位置的基因?qū)儆谀硞€(gè)指定通路巩剖。綠色有波動(dòng)的曲線表示富集分?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è)的話說就是先生成一個(gè)零假設(shè)的數(shù)據(jù)分布，然后觀察實(shí)際數(shù)據(jù)在這個(gè)零假設(shè)分布下稳衬，是不是在尾端霞捡。
一些問題
統(tǒng)計(jì)檢驗(yàn)的能力是有限的，所以還有很多問題存在解決薄疚。
- 我們希望找到生物學(xué)顯著的基因碧信，但是生物學(xué)顯著和統(tǒng)計(jì)顯著兩者并不是完全相關(guān)
- 無論是ORA還是FCS都對(duì)背景（也就是這個(gè)物種一共有多少基因）有要求赊琳，但是隨著我們的研究深入，基因數(shù)量會(huì)改變音婶。有些軟件會(huì)直接設(shè)置一個(gè)很大的背景數(shù)慨畸，從而讓p值很顯著，然后我們就開心地用他們的結(jié)果衣式。
- 有些基因沒有注釋寸士，也就是注釋缺失，處理方法就是扔（歡迎拍磚）碴卧。
- 有一些注釋項(xiàng)是其他項(xiàng)的子集弱卡。
實(shí)戰(zhàn)：用clusterProfiler做富集分析
clusterProfiler是Y叔良心之作（當(dāng)然他的良心之作還有很多），目前支持KEGG在線拉數(shù)據(jù)住册，支持DAVID婶博，支持Broad iNSTITUTE Molecular Signatures Databases， 支持GSEA荧飞。
安裝：
source("https://bioconductor.org/biocLite.R")
biocLite('clusterProfiler')
這里簡單舉例如何使用凡人，更多內(nèi)容見軟件文檔，當(dāng)然這個(gè)部分你可以直接過掉叹阔，因?yàn)槲抑皇歉鳼叔的代碼敲了一遍而已挠轴。
加載差異表達(dá)基因，我這里偷懶就隨機(jī)挑一些基因名（來自于DOSE包耳幢，Disease Ontology Semantic and Enrichment analysis ）出來了岸晦。
library(org.Hs.eg.db)
data(geneList)
gene <- names(geneList)[abs(geneList) > 2]
head(gene)
[1] "4312" "8318" "10874" "55143" "55388" "991"
Y叔為了展示他能夠處理不同命名方式的ID，用bitr（來自于clusterProfiler)進(jìn)行生物學(xué)ID轉(zhuǎn)換
gene.df <- bitr(gene, fromType = "ENTREZID",
toType = c("ENSEMBL", "SYMBOL"),
OrgDb = org.Hs.eg.db)
head(gene.df)
ENTREZID ENSEMBL SYMBOL
1 4312 ENSG00000196611 MMP1
2 8318 ENSG00000093009 CDC45
3 10874 ENSG00000109255 NMU
4 55143 ENSG00000134690 CDCA8
5 55388 ENSG00000065328 MCM10
6 991 ENSG00000117399 CDC20
clusterProfiler: ORA
然后對(duì)不同命名的ID都做ORA富集分析睛藻。
ego <- enrichGO(gene = gene,
universe = names(geneList),
OrgDb = org.Hs.eg.db,
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05)
ego2 <- enrichGO(gene = gene.df$ENSEMBL,
OrgDb = org.Hs.eg.db,
keytype = 'ENSEMBL',
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05)
ego3 <- enrichGO(gene = gene.df$SYMBOL,
OrgDb = org.Hs.eg.db,
keytype = 'SYMBOL',
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05)
結(jié)果會(huì)有如下
ID： 基因本體論的ID
Description： 描述
GeneRatio： 在GO詞條所占的比例
BgRatio ：在背景所占的比例
pvalue： 假設(shè)是正確但是被拒絕的概率
p.adjust： 采用BH方法進(jìn)行多重試驗(yàn)p值校正
qvalue： Q值=被拒絕但卻是正確的概率
geneID： 列出在這個(gè)GO詞條下的我們提供的基因
count： 數(shù)量启上，如果沒有約分，就是GeneRatio的分子了店印。
所以從GeneRatio：24/199冈在， BgRatio：231/11711可以反推出下表：
| 類別 | gene list | Genome |
|---|---|---|
| in anno group | 24 | 199 |
| not in anno group | 231 | 11711 |
| 255 | 11910 |
最后再畫一張美美的點(diǎn)圖,根據(jù)GeneRatio所做。
dotplot(ego, showCategory=30)
clusterProfiler: GSEC
clusterProfiler支持GSEC按摘，而且很好用讥邻。
gsecc <- gseGO(geneList=geneList, ont="CC", OrgDb=org.Hs.eg.db, verbose=F)
head(summary(gsecc)
gseaplot(gsecc, geneSetID="GO:0000779"
