GO、GSEA富集分析一網(wǎng)打進(jìn)
GO、KEGG、GSEA富集分析，已經(jīng)注釋好的數(shù)據(jù)，自己注釋的數(shù)據(jù)，都在這里
NGS系列文章包括NGS基礎(chǔ)、轉(zhuǎn)錄組分析?（Nature重磅綜述|關(guān)于RNA-seq你想知道的全在這）、ChIP-seq分析?（ChIP-seq基本分析流程）、單細(xì)胞測(cè)序分析?(重磅綜述：三萬(wàn)字長(zhǎng)文讀懂單細(xì)胞RNA測(cè)序分析的最佳實(shí)踐教程 （原理、代碼和評(píng)述）)、DNA甲基化分析、重測(cè)序分析、GEO數(shù)據(jù)挖掘（典型醫(yī)學(xué)設(shè)計(jì)實(shí)驗(yàn)GEO數(shù)據(jù)分析 (step-by-step) - Limma差異分析、火山圖、功能富集）等內(nèi)容。
富集分析是生物信息分析中快速了解目標(biāo)基因或目標(biāo)區(qū)域功能傾向性的最重要方法之一。其中代表性的計(jì)算方式有兩種：
一是基于篩選的差異基因，采用超幾何檢驗(yàn)判斷上調(diào)或下調(diào)基因在哪些GO或KEGG或其它定義的通路富集。假設(shè)背景基因數(shù)目為m，背景基因中某一通路pathway中注釋的基因有n個(gè)；上調(diào)基因有k個(gè)，上調(diào)基因中落于通路pathway的數(shù)目為l。簡(jiǎn)單來(lái)講就是比較l/k是否顯著高于n/m，即上調(diào)基因中落在通路pathway的比例是否高于背景基因在這一通路的比例。(實(shí)際計(jì)算時(shí)，是算的odds ratio的差異，l/(k-l) vs (n-l)/(m-k-n+l))。這就是常說(shuō)的GO富集分析或KEGG富集分析，可以做的工具很多，GOEAST是其中一個(gè)最好用的在線功能富集分析工具，數(shù)據(jù)庫(kù)更新實(shí)時(shí)，操作簡(jiǎn)單，并且可以直接用之前介紹的方法繪制DotPlot。
另一種方式是不硬篩選差異基因，而是對(duì)其根據(jù)表達(dá)量或與表型的相關(guān)度排序，然后判斷對(duì)應(yīng)的基因集是否傾向于落在有序列表的頂部或底部，從而判斷基因集合對(duì)表型差異的影響和篩選有影響的基因子集。這叫GSEA富集分析，注釋信息可以是GO，KEGG，也可以是其它任何符合格式的信息。GSEA富集分析 - 界面操作詳細(xì)講述了GSEA分析的原理、可視化操作和結(jié)果解讀 （一文掌握GSEA，超詳細(xì)教程）。
R包clusterProfiler囊括了前面提到的兩種富集分析方法。并且不只支持GO、KEGG數(shù)據(jù)庫(kù)，還支持Disease Ontology、MsEH enrichment analysis、Reactome通路分析等。具體可見(jiàn)軟件的文檔頁(yè) https://guangchuangyu.github.io/clusterProfiler/。
這么強(qiáng)大的工具，學(xué)習(xí)起來(lái)的路子卻不是一帆風(fēng)順，最開(kāi)始的攔路虎是軟件的安裝，系統(tǒng)較老配合上軟件包更新較快(作者勤快也是問(wèn)題),導(dǎo)致經(jīng)常安裝的是舊版本，用起來(lái)會(huì)遇到不少坑。直到有了conda，安裝再也不是問(wèn)題。解決了動(dòng)態(tài)庫(kù)依賴(lài)后，可以在Github安裝最新的開(kāi)發(fā)版本。
另外一個(gè)是文檔較少，在R終端，直接使用help命令查看到的使用提示信息較少，寥寥幾句，看過(guò)總覺(jué)得不踏實(shí)。~在線文檔頁(yè)內(nèi)容少、更新慢。這是最開(kāi)始學(xué)習(xí)時(shí)遇到的問(wèn)題，這次秉著負(fù)責(zé)的精神，又重新讀了文檔頁(yè)，發(fā)覺(jué)不需要再寫(xiě)一遍了，內(nèi)容挺全的，但有幾個(gè)地方需要更新下。自己對(duì)著文檔頁(yè)核對(duì)了下之前寫(xiě)的程序，再補(bǔ)充幾點(diǎn)。
GO富集分析
首先還是列一個(gè)完整的例子。輸入最好是用ENTREZ ID，值比較固定，不建議使用GeneSymbol，容易匹配出問(wèn)題 （Excel改變了你的基因名，30% 相關(guān)Nature文章受影響，NCBI也受波及）。
entrezID_text <- "4312
8318
10874
55143
55388
991
6280
2305
9493
1062
3868
4605
9833
9133
6279
10403
8685
597
7153
23397"
entrezID <- read.table(text=entrezID_text, header=F)
head(entrezID)V1 ? ? ? ? ? ?
1 ? ?4312 ? ? ? ? ? ?
2 ? ?8318 ? ? ? ? ? ?
3 ? ?10874 ? ? ? ? ? ?
4 ? ?55143 ? ? ? ? ? ?
5 ? ?55388 ? ? ? ? ? ?
6 ? ?991轉(zhuǎn)換為向量
entrezID <- entrezID$V1
head(entrezID)[1] 4312 8318 10874 55143 55388 991# 這里的ENTREZ ID是從clusterProfiler里面提取是，是人的，
# 所以用了人的注釋庫(kù), org.Hs.eg.db# 這個(gè)庫(kù)半年左右更新一次，也有可能漏更新，還是自己準(zhǔn)備數(shù)據(jù)的靠譜
library(org.Hs.eg.db)開(kāi)始富集分析
# readable=T: 原文檔無(wú)這個(gè)參數(shù)，使用的是setReadble函數(shù)
MF <- enrichGO(entrezID, "org.Hs.eg.db", ont = "MF", keytype = "ENTREZID",
? ? ? ? ?pAdjustMethod = "BH", pvalueCutoff ?= 0.05,
? ? ? ? ?qvalueCutoff ?= 0.1, readable=T)
head(summary(MF))
去除冗余度高的條目?
(參考http://guangchuangyu.github.io/2015/10/use-simplify-to-remove-redundancy-of-enriched-go-terms/)
result <- simplify(MF, cutoff=0.7, by="p.adjust", select_fun=min)
# 去除前
dim(MF)
# [1] 367 9
# 去除后
dim(result)
[1] 142 9繪制泡泡圖 (這里還有一個(gè)可以直接在線繪制泡泡圖的網(wǎng)站：http://www.ehbio.com/ImageGP/)
dotplot(result, showCategory = 10)繪制網(wǎng)絡(luò)圖 (邊的寬度代表兩個(gè)富集的Term共有的基因數(shù)目，點(diǎn)大小代表?xiàng)l目?jī)?nèi)基因數(shù)目多少，顏色代表P-value，值越小越紅；如果想改變網(wǎng)絡(luò)布局，參考igraph文檔)
enrichMap(result, vertex.label.cex=1.2, layout=igraph::layout.kamada.kawai)另外一種網(wǎng)絡(luò)圖由cnetplot函數(shù)獲得，可以映射基因的表達(dá)量。
# geneList為一個(gè)vector，每個(gè)元素的名字為基因名，值為FoldChange，用于可視化點(diǎn)。
# > str(geneList)
# Named num [1:12495] 4.57 4.51 4.42 4.14 3.88 ... - attr(*, "names")= chr [1:12495] "4312" "8318" "10874" "55143" ...
# 這個(gè)geneList怎么獲得的，會(huì)在后面的GSEA分析時(shí)提到
cnetplot(result, foldChange=geneList)自己嘗試了下，展示的有些亂，需要調(diào)整字體和顯示的條目多少。故盜圖展示如下便于解釋?zhuān)蚺c其被注釋的條目連線，點(diǎn)的顏色代表表達(dá)變化，圈的大小代表對(duì)應(yīng)注釋內(nèi)基因數(shù)目多少。
如果想自己調(diào)整圖的布局，還是建議把輸出結(jié)果轉(zhuǎn)換為Cytoscape(點(diǎn)擊查看視頻教程)可以識(shí)別的兩列表格形式(如下)，再賦值不同的屬性就可以了。
awk 'BEGIN{OFS=FS="\t"}{if(FNR==1) print "Gene\tTerm"; else {split($8,geneL,"/"); for(i in geneL) print geneL[i],$2}}' MF | head
Gene ? ?Term
PRDX6 ? ?cell adhesion molecule binding
MFGE8 ? ?cell adhesion molecule binding
FSCN1 ? ?cell adhesion molecule binding
ATXN2L ? ?cell adhesion molecule binding
YWHAZ ? ?cell adhesion molecule binding
CTNNA2 ? ?cell adhesion molecule binding
ADAM15 ? ?cell adhesion molecule binding
LDHA ? ?cell adhesion molecule binding
PKM ? ?cell adhesion molecule bindingKEGG富集分析
輸入Gene ID的格式和類(lèi)型與enrichGO一致。參考
http://guangchuangyu.github.io/2015/02/kegg-enrichment-analysis-with-latest-online-data-using-clusterprofiler/。
# clusterProfiler3.4.4版本是沒(méi)有readable參數(shù)的，原文檔有
kk <- enrichKEGG(entrezID, organism="hsa", pvalueCutoff=0.05, pAdjustMethod="BH",
? ? ? ? ? ? ? ? qvalueCutoff=0.1)輸出結(jié)果的格式和可視化方式與上面GO富集一致，不再贅述。
另外一個(gè)沒(méi)有解決的問(wèn)題是setReadable函數(shù)的使用 (用測(cè)試文檔提供的數(shù)據(jù)集報(bào)出如下錯(cuò)誤)
setReadable(kk, "org.Hs.eg.db")
Error in EXTID2NAME(OrgDb, genes, keytype) : keytype is not supported...
# 即便是如下操作也沒(méi)有作用
a = bitr_kegg(names(geneList), fromType = 'ncbi-geneid', toType="kegg", organism="hsa")
# Warning message:
# In bitr_kegg(names(geneList), fromType = "ncbi-geneid", toType = "kegg", ?:
# ?0.77% of input gene IDs are fail to map...
kk <- enrichKEGG(a$kegg, organism="hsa", keyType="kegg", pvalueCutoff=0.05, pAdjustMethod="BH",
? ? ? ? ? ? ? ? ?qvalueCutoff=0.1)
setReadable(kk, org.Hs.eg.db, key)
# ?Error in EXTID2NAME(OrgDb, genes, keytype) : keytype is not supported...經(jīng)過(guò)多次嘗試發(fā)現(xiàn)，可以這么解決
setReadable(kk, ?org.Hs.eg.db, keytype="ENTREZID")為什么會(huì)有這個(gè)問(wèn)題呢？setReadable中自動(dòng)判斷keytype的語(yǔ)句是
if (keytype == "auto") {
? ?keytype <- x@keytype
? ?if (keytype == "UNKNOWN") {
? ? ? ?stop("can't determine keytype automatically; need to set 'keytype' explicitly...")
? ? ? ?}
}根據(jù)富集結(jié)果中定義的keytype，也就是enrichKEGG函數(shù)中設(shè)定的keyType的值來(lái)定的。而setReadable不支持默認(rèn)的keyType=kegg。
沒(méi)有測(cè)試小鼠，可能需要設(shè)置不同的keytype值。
另外對(duì)擬南芥來(lái)說(shuō)，分析之前需要先把Entrez ID轉(zhuǎn)換為kegg再用上述命令做富集分析
entrezID <- bitr_kegg(entrezID, fromType='ncbi-geneid', toType='kegg', organism="ath")
kk <- enrichKEGG(entrezID$kegg, organism="ath", pvalueCutoff=0.05, pAdjustMethod="BH",
? ? ? ? ? ? ? ? qvalueCutoff=0.1)
# 這個(gè)setreadble是可以轉(zhuǎn)換成功的
result <- setReadable(kk, "org.At.tair.db", keytype="TAIR")GSEA分析
GSEA的解釋和介紹見(jiàn)GSEA富集分析 - 界面操作。
注意讀入的基因列表是要按照表達(dá)差異降序排列 (升序也可以，相當(dāng)于樣品做了對(duì)調(diào))。這里排序方式可以是表達(dá)差異，也可以是其它方式，只要方便解釋即可，即從上到下，或從前到后，基因?qū)Ρ硇偷呢暙I(xiàn)有一致的變化趨勢(shì)就好。不同的排序參數(shù)和排序方式需要不同的對(duì)結(jié)果的解釋。
id_with_fc = "ID;FC
4312;2
8318;3
10874;4
55143;5
55388;6
991;7"
id_with_fc <- read.table(text=id_with_fc, header=T, sep=";")
id_with_fc2 <- id_with_fc[,2]
names(id_with_fc2) <- id_with_fc[,1]
# 排序是必須的，記住排序方式
id_with_fc2 <- sort(id_with_fc2, decreasing=T)
gsecc <- gseGO(geneList=id_with_fc2, ont="CC", OrgDb=org.Hs.eg.db, verbose=F)
# 昨天測(cè)試了其它數(shù)據(jù)，參數(shù)無(wú)問(wèn)題。這里沒(méi)有實(shí)際運(yùn)行，盜用數(shù)據(jù)，每列的解釋見(jiàn)本段開(kāi)頭的文章
head(as.data.frame(gsecc))
## ? ? ? ? ? ? ? ? ? ?ID ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?Description setSize
## GO:0031982 GO:0031982 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?vesicle ? ?2880
## GO:0031988 GO:0031988 ? ? ? ? ? ? ? ? membrane-bounded vesicle ? ?2791
## GO:0005576 GO:0005576 ? ? ? ? ? ? ? ? ? ? extracellular region ? ?3296
## GO:0065010 GO:0065010 extracellular membrane-bounded organelle ? ?2220
## GO:0070062 GO:0070062 ? ? ? ? ? ? ? ? ? ?extracellular exosome ? ?2220
## GO:0044421 GO:0044421 ? ? ? ? ? ? ? ?extracellular region part ? ?2941
## ? ? ? ? ? ?enrichmentScore ? ? ? NES ? ? ?pvalue ? p.adjust ? ?qvalues
## GO:0031982 ? ? ?-0.2561837 -1.222689 0.001002004 0.03721229 0.02816364
## GO:0031988 ? ? ?-0.2572169 -1.226003 0.001007049 0.03721229 0.02816364
## GO:0005576 ? ? ?-0.2746489 -1.312485 0.001009082 0.03721229 0.02816364
## GO:0065010 ? ? ?-0.2570342 -1.222048 0.001013171 0.03721229 0.02816364
## GO:0070062 ? ? ?-0.2570342 -1.222048 0.001013171 0.03721229 0.02816364
## GO:0044421 ? ? ?-0.2744658 -1.310299 0.001014199 0.03721229 0.02816364
# 繪制GSEA圖
gseaplot(gsecc, geneSetID="GO:0000779")自定義數(shù)據(jù)集分析
如果想用clusterProfiler的函數(shù)對(duì)自己注釋的數(shù)據(jù)進(jìn)行功能富集分析或GSEA分析，需要提供如下格式的注釋數(shù)據(jù)。后續(xù)分析就類(lèi)似了。
self_anno <- "ont;gene
KEGG_GLYCOLYSIS_GLUCONEOGENESIS;gene1
KEGG_GLYCOLYSIS_GLUCONEOGENESIS;gene2
KEGG_GLYCOLYSIS_GLUCONEOGENESIS;gene3
KEGG_GLYCOLYSIS;gene1
KEGG_GLYCOLYSIS;gene4
KEGG_CYP;gene5"
self_anno <- read.table(text=self_anno, header=T, sep=";", quote="")
# 沒(méi)具體看代碼怎么寫(xiě)的，保險(xiǎn)期間，設(shè)置跟示例一樣的列名字
colnames(self_anno) <- c("ont", "gene")
geneL <- c("gene1", "gene2", "gene4")
# self_enrich與之前enrichGO的輸出結(jié)果格式一致
self_enrich <- enricher(geneL, TERM2GENE=self_anno)
# self_gsea與之前gseGO的輸出結(jié)果格式一致
self_gsea <- GSEA(geneL, TERM2GENE=self_anno, verbose=F)Reference
https://guangchuangyu.github.io/clusterProfiler/
http://guangchuangyu.github.io/2016/01/go-analysis-using-clusterprofiler/
http://igraph.org/c/doc/igraph-Layout.html
你可能還想看
往期精品(點(diǎn)擊圖片直達(dá)文字對(duì)應(yīng)教程)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
后臺(tái)回復(fù)“生信寶典福利第一波”或點(diǎn)擊閱讀原文獲取教程合集
