亚洲日本中文字幕乱码在线 ,国产成人三级在线视频,久久日大香蕉,国产激情高清视频,97精品自拍,91作爱视频,水密桃网站,人人操超碰在线观看

GO、KEGG、GSEA富集分析，已經(jīng)注釋好的數(shù)據(jù)，自己注釋的數(shù)據(jù)，都在這里

NGS系列文章包括NGS基礎(chǔ)、轉(zhuǎn)錄組分析?（Nature重磅綜述|關(guān)于RNA-seq你想知道的全在這）、ChIP-seq分析?（ChIP-seq基本分析流程）、單細胞測序分析?(重磅綜述：三萬字長文讀懂單細胞RNA測序分析的最佳實踐教程（原理、代碼和評述）)、DNA甲基化分析、重測序分析、GEO數(shù)據(jù)挖掘（典型醫(yī)學設(shè)計實驗GEO數(shù)據(jù)分析 (step-by-step) - Limma差異分析、火山圖、功能富集）等內(nèi)容。

富集分析是生物信息分析中快速了解目標基因或目標區(qū)域功能傾向性的最重要方法之一。其中代表性的計算方式有兩種：

一是基于篩選的差異基因，采用超幾何檢驗判斷上調(diào)或下調(diào)基因在哪些GO或KEGG或其它定義的通路富集。假設(shè)背景基因數(shù)目為m，背景基因中某一通路pathway中注釋的基因有n個；上調(diào)基因有k個，上調(diào)基因中落于通路pathway的數(shù)目為l。簡單來講就是比較l/k是否顯著高于n/m，即上調(diào)基因中落在通路pathway的比例是否高于背景基因在這一通路的比例。(實際計算時，是算的odds ratio的差異，l/(k-l) vs (n-l)/(m-k-n+l))。這就是常說的GO富集分析或KEGG富集分析，可以做的工具很多，GOEAST是其中一個最好用的在線功能富集分析工具，數(shù)據(jù)庫更新實時，操作簡單，并且可以直接用之前介紹的方法繪制DotPlot。

另一種方式是不硬篩選差異基因，而是對其根據(jù)表達量或與表型的相關(guān)度排序，然后判斷對應(yīng)的基因集是否傾向于落在有序列表的頂部或底部，從而判斷基因集合對表型差異的影響和篩選有影響的基因子集。這叫GSEA富集分析，注釋信息可以是GO，KEGG，也可以是其它任何符合格式的信息。GSEA富集分析 - 界面操作詳細講述了GSEA分析的原理、可視化操作和結(jié)果解讀（一文掌握GSEA，超詳細教程）。

R包clusterProfiler囊括了前面提到的兩種富集分析方法。并且不只支持GO、KEGG數(shù)據(jù)庫，還支持Disease Ontology、MsEH enrichment analysis、Reactome通路分析等。具體可見軟件的文檔頁 https://guangchuangyu.github.io/clusterProfiler/。

這么強大的工具，學習起來的路子卻不是一帆風順，最開始的攔路虎是軟件的安裝，系統(tǒng)較老配合上軟件包更新較快(作者勤快也是問題),導致經(jīng)常安裝的是舊版本，用起來會遇到不少坑。直到有了conda，安裝再也不是問題。解決了動態(tài)庫依賴后，可以在Github安裝最新的開發(fā)版本。

另外一個是文檔較少，在R終端，直接使用help命令查看到的使用提示信息較少，寥寥幾句，看過總覺得不踏實。~在線文檔頁內(nèi)容少、更新慢。這是最開始學習時遇到的問題，這次秉著負責的精神，又重新讀了文檔頁，發(fā)覺不需要再寫一遍了，內(nèi)容挺全的，但有幾個地方需要更新下。自己對著文檔頁核對了下之前寫的程序，再補充幾點。

GO富集分析

首先還是列一個完整的例子。輸入最好是用ENTREZ ID，值比較固定，不建議使用GeneSymbol，容易匹配出問題（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里面提取是，是人的，

# 所以用了人的注釋庫, org.Hs.eg.db# 這個庫半年左右更新一次，也有可能漏更新，還是自己準備數(shù)據(jù)的靠譜

library(org.Hs.eg.db)

開始富集分析

# readable=T: 原文檔無這個參數(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

繪制泡泡圖 (這里還有一個可以直接在線繪制泡泡圖的網(wǎng)站：http://www.ehbio.com/ImageGP/)

dotplot(result, showCategory = 10)

繪制網(wǎng)絡(luò)圖 (邊的寬度代表兩個富集的Term共有的基因數(shù)目，點大小代表條目內(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ù)獲得，可以映射基因的表達量。

# 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)

自己嘗試了下，展示的有些亂，需要調(diào)整字體和顯示的條目多少。故盜圖展示如下便于解釋，基因與其被注釋的條目連線，點的顏色代表表達變化，圈的大小代表對應(yīng)注釋內(nèi)基因數(shù)目多少。

如果想自己調(diào)整圖的布局，還是建議把輸出結(jié)果轉(zhuǎn)換為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 binding

KEGG富集分析

輸入Gene ID的格式和類型與enrichGO一致。參考

http://guangchuangyu.github.io/2015/02/kegg-enrichment-analysis-with-latest-online-data-using-clusterprofiler/。

# clusterProfiler3.4.4版本是沒有readable參數(shù)的，原文檔有
kk <- enrichKEGG(entrezID, organism="hsa", pvalueCutoff=0.05, pAdjustMethod="BH",
 ? ? ? ? ? ? ? ? qvalueCutoff=0.1)

輸出結(jié)果的格式和可視化方式與上面GO富集一致，不再贅述。

另外一個沒有解決的問題是setReadable函數(shù)的使用 (用測試文檔提供的數(shù)據(jù)集報出如下錯誤)

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...

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)過多次嘗試發(fā)現(xiàn)，可以這么解決

setReadable(kk, ?org.Hs.eg.db, keytype="ENTREZID")

為什么會有這個問題呢？setReadable中自動判斷keytype的語句是

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的值來定的。而setReadable不支持默認的keyType=kegg。

沒有測試小鼠，可能需要設(shè)置不同的keytype值。

另外對擬南芥來說，分析之前需要先把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)
# 這個setreadble是可以轉(zhuǎn)換成功的
result <- setReadable(kk, "org.At.tair.db", keytype="TAIR")

GSEA分析

GSEA的解釋和介紹見GSEA富集分析 - 界面操作。

注意讀入的基因列表是要按照表達差異降序排列 (升序也可以，相當于樣品做了對調(diào))。這里排序方式可以是表達差異，也可以是其它方式，只要方便解釋即可，即從上到下，或從前到后，基因?qū)Ρ硇偷呢暙I有一致的變化趨勢就好。不同的排序參數(shù)和排序方式需要不同的對結(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)

# 昨天測試了其它數(shù)據(jù)，參數(shù)無問題。這里沒有實際運行，盜用數(shù)據(jù)，每列的解釋見本段開頭的文章
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ù)對自己注釋的數(shù)據(jù)進行功能富集分析或GSEA分析，需要提供如下格式的注釋數(shù)據(jù)。后續(xù)分析就類似了。

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="")

# 沒具體看代碼怎么寫的，保險期間，設(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)