富集分析一网打进
歡迎關注微信公眾號生信寶典：http://mp.weixin.qq.com/s/d1KCETQZ88yaOLGwAtpWYg
富集分析是生物信息分析中快速了解目標基因或目標區域功能傾向性的最重要方法之一。其中代表性的計算方式有兩種：
一是基于篩選的差異基因，采用超幾何檢驗判斷上調或下調基因在哪些GO或KEGG或其它定義的通路富集。假設背景基因數目為m，背景基因中某一通路pathway中注釋的基因有n個；上調基因有k個，上調基因中落于通路pathway的數目為l。簡單來講就是比較l/k是否顯著高于n/m，即上調基因中落在通路pathway的比例是否高于背景基因在這一通路的比例。(實際計算時，是算的odds ratio的差異，l/(k-l) vs (n-l)/(m-k-n+l))。這就是常說的GO富集分析或KEGG富集分析，可以做的工具很多，GOEAST是其中一個最好用的在線功能富集分析工具，數據庫更新實時，操作簡單，并且可以直接用之前介紹的方法繪制DotPlot。
另一種方式是不硬篩選差異基因，而是對其根據表達量或與表型的相關度排序，然后判斷對應的基因集是否傾向于落在有序列表的頂部或底部，從而判斷基因集合對表型差異的影響和篩選有影響的基因子集。這叫GSEA富集分析，注釋信息可以是GO，KEGG，也可以是其它任何符合格式的信息。GSEA富集分析 - 界面操作詳細講述了GSEA分析的原理、可視化操作和結果解讀。
微信公眾號biobabble的主創Y叔寫有9個Bioconductor包，其中兩個DOSE和clusterProfiler囊括了前面提到的兩種富集分析方法。并且不只支持GO、KEGG數據庫，還支持Disease Ontology、MsEH enrichment analysis、Reactome通路分析等。具體可見其公眾號，或軟件的文檔頁 https://guangchuangyu.github.io/clusterProfiler/。
這么強大的工具，學習起來的路子卻不是一帆風順，最開始的攔路虎是軟件的安裝，系統較老配合上軟件包更新較快，導致經常安裝的是舊版本，用起來會遇到不少坑。直到有了conda，安裝再也不是問題。解決了動態庫依賴后，可以在Github安裝最新的修改。
另外一個是文檔較少，在R終端，直接使用help命令查看到的使用提示信息較少，寥寥幾句，看過總覺得不踏實。~在線文檔頁內容少、更新慢~。這是最開始學習時遇到的問題，這次秉著負責的精神，又重新讀了文檔頁，發覺不需要再寫一遍了，內容挺全的，主要是這一頁http://guangchuangyu.github.io/2016/01/go-analysis-using-clusterprofiler/)，但有幾個地方需要更新下。自己對著文檔頁核對了下之前寫的程序，再補充幾點。
GO富集分析
首先還是列一個完整的例子。輸入最好是用ENTREZ ID，值比較固定，不建議使用GeneSymbol，容易匹配出問題。
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轉換為向量
entrezID <- entrezID$V1 head(entrezID) [1] 4312 8318 10874 55143 55388 991 # 這里的ENTREZ ID是從clusterProfiler里面提取是，是人的， # 所以用了人的注釋庫, org.Hs.eg.db library(org.Hs.eg.db)開始富集分析
# readable=T: 原文檔無這個參數，使用的是setReadble函數 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繪制泡泡圖 （DotPlot Dotplot2）
dotplot(result, showCategory = 10)繪制網絡圖(邊的寬度代表兩個富集的Term共有的基因數目，點大小代表條目內基因數目多少，顏色代表P-value，值越小越紅；如果想改變網絡布局，參考igraph文檔)
enrichMap(result, vertex.label.cex=1.2, layout=igraph::layout.kamada.kawai)另外一種網絡圖由cnetplot函數獲得，可以映射基因的表達量。
# geneList為一個vector，每個元素的名字為基因名，值為FoldChange，用于可視化點。 # > 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" ...# 這個geneList怎么獲得的，會在后面的GSEA分析時提到 cnetplot(result, foldChange=geneList)自己嘗試了下，展示的有些亂，需要調整字體和顯示的條目多少。故盜圖展示如下便于解釋，基因與其被注釋的條目連線，點的顏色代表表達變化，圈的大小代表對應注釋內基因數目多少。
如果想自己調整圖的布局，還是建議把輸出結果轉換為Cytoscape可以識別的兩列表格形式(如下)，再賦值不同的屬性就可以了。
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的格式和類型與enrichGO一致。參考 http://guangchuangyu.github.io/2015/02/kegg-enrichment-analysis-with-latest-online-data-using-clusterprofiler/。
# clusterProfiler3.4.4版本是沒有readable參數的，原文檔有 kk <- enrichKEGG(entrezID, organism="hsa", pvalueCutoff=0.05, pAdjustMethod="BH", qvalueCutoff=0.1)輸出結果的格式和可視化方式與上面GO富集一致，不再贅述。
另外一個沒有解決的問題是setReadable函數的使用 (用測試文檔提供的數據集報出如下錯誤)
setReadable(kk, "org.Hs.eg.db") Error in EXTID2NAME(OrgDb, genes, keytype) : keytype is not supported...# 即便是如下操作也沒有作用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...enrichKEGG(a$kegg, organism="hsa", keyType="kegg", pvalueCutoff=0.05, pAdjustMethod="BH", qvalueCutoff=0.1)setReadable(kk, org.Hs.eg.db)# Error in EXTID2NAME(OrgDb, genes, keytype) : keytype is not supported...經過多次嘗試發現，可以這么解決
setReadable(kk, org.Hs.eg.db, keytype="ENTREZID")為什么會有這個問題呢？setReadable中自動判斷keytype的語句是
if (keytype == "auto") {keytype <- x@keytypeif (keytype == "UNKNOWN") {stop("can't determine keytype automatically; need to set 'keytype' explicitly...")} }根據富集結果中定義的keytype，也就是enrichKEGG函數中設定的keyType的值來定的。而setReadable不支持默認的keyType=kegg。
沒有測試小鼠，可能需要設置不同的keytype值。
對擬南芥來說，分析之前需要先把Entrez ID轉換為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) # 這個setreadble是可以轉換成功的，這里的keytype是TAIR result <- setReadable(kk, "org.At.tair.db", keytype="TAIR")GSEA分析
GSEA的解釋和介紹見GSEA富集分析 - 界面操作。
注意讀入的基因列表是要按照表達差異降序排列 (升序也可以，相當于樣品做了對調)。這里排序方式可以是表達差異，也可以是其它方式，只要方便解釋即可，即從上到下，或從前到后，基因對表型的貢獻有一致的變化趨勢就好。不同的排序參數和排序方式需要不同的對結果的解釋。
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)# 昨天測試了其它數據，參數無問題。這里沒有實際運行，盜用數據，每列的解釋見本段開頭的文章 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")自定義數據集分析
如果想用clusterProfiler的函數對自己注釋的數據進行功能富集分析或GSEA分析，需要提供如下格式的注釋數據。后續分析就類似了。
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="")# 沒具體看代碼怎么寫的，保險期間，設置跟示例一樣的列名字 colnames(self_anno) <- c("ont", "gene")geneL <- c("gene1", "gene2", "gene4")# self_enrich與之前enrichGO的輸出結果格式一致 self_enrich <- enricher(geneL, TERM2GENE=self_anno)# self_gsea與之前gseGO的輸出結果格式一致 self_gsea <- GSEA(geneL, TERM2GENE=self_anno, verbose=F)九天學會轉錄組高級分析
經過緊張的籌備，生信寶典團隊要開培訓課了，第一期是大家最為關注的轉錄組分析，實行三段式培訓，集中講解實戰(2天)-自行練習(5天)-再講解答疑考核(2天)。點擊原文可查看具體信息。
歡迎大家咨詢、報名和提出建議。
Reference
生信寶典，生物信息學習系列教程，轉錄組，宏基因組，外顯子組，R作圖，Python學習，Cytoscape視頻教程
http://mp.weixin.qq.com/s/d1KCETQZ88yaOLGwAtpWYg
生信寶典，最好的生物信息培訓課程，培訓課程資料
www.ehbio.com/Training
總結
- 上一篇： ON1 photo raw 2021（p
- 下一篇： 施一公：如何提高英文的科研写作能力