https://zhuanlan.zhihu.com/p/597338272
https://zhuanlan.zhihu.com/p/569459288
http://www.reibang.com/p/25386890a751
http://www.reibang.com/p/0b80b24b0d03
http://www.bio-info-trainee.com/6846.html
https://blog.csdn.net/qq_42090739/article/details/123246849
目錄
收起
1 用 DOSE 包里的基因list為例
2 clusterProfiler包進行GO富集分析
2.1 加載包
2.2 enrichGO函數(shù)富集分析
3 enrichplot包繪圖
3.1 barplot
3.2 dotplot
4 ggplot2包繪圖
4.1 計算Enrichment Factor
4.2 BP\MF\CC各取排名前10的term
4.3 ggplot2畫圖
4.4 ggplot2美化
4.5 ggplot2分屏繪圖
4.5 ggplot2調(diào)整排序
5 GOBubble圖
5.1 準備circle_dat數(shù)據(jù)
5.2 繪制GOBubble圖
6 GOCircle圖
7 GOHeat熱圖
8 GOChord弦圖
9 樹狀圖
通過某些方法（如差異基因、相關(guān)基因等）獲得一批基因后昧互，可以進行功能富集分析（包括GO和KEGG富集分析）航瞭，以了解這些基因的主要功能和涉及信號通路有哪些袋倔。
本文用clusterProfiler包GO富集分析，分別用enrichplot和ggplot2繪圖诅岩，并用GOplot包繪制Bubble圖讳苦、Circle圖、Heat圖吩谦、Chord圖鸳谜。
1 用DOSE包里的基因list為例
library("DOSE")
data(geneList, package = "DOSE")
g_list <- names(geneList)[1:100]
head(g_list)
[1] "4312" "8318" "10874" "55143" "55388" "991"
從DOSE包里取出了100個基因的EntrezID。
2 clusterProfiler包進行GO富集分析
2.1 加載包
library("clusterProfiler")
library("org.Hs.eg.db")
2.2 enrichGO函數(shù)富集分析
eG <- enrichGO(gene = g_list, #需要分析的基因的EntrezID
OrgDb = org.Hs.eg.db, #人基因數(shù)據(jù)庫
pvalueCutoff =0.01, #設(shè)置pvalue界值
qvalueCutoff = 0.01, #設(shè)置qvalue界值(FDR校正后的p值）
ont="all", #選擇功能富集的類型式廷，可選BP咐扭、MF、CC滑废，這里選擇all蝗肪。
readable =T)
富集分析的類型可選BP（biological process）、MF（molecular function）策严、CC（Cellular Component）穗慕。
輸出結(jié)果為txt文件
write.table(eG,file="eG.txt", sep="\t", quote=F, row.names = F)
結(jié)果如下：
[圖片上傳失敗...(image-c1d626-1679977340335)]
<figcaption style="color: rgb(153, 153, 153); font-size: 0.9em; line-height: 1.5; margin-top: calc(0.666667em); padding: 0px 1em; text-align: center;">GO</figcaption>
3 enrichplot包繪圖
enrichplot包可以直接用上文enrichGO分析出的eG文件繪圖。
3.1 barplot
pdf(file="eGO_barplot.pdf",width = 8,height = 10)
barplot(eG, x = "GeneRatio", color = "p.adjust", #默認參數(shù)（x和color可以根據(jù)eG里面的內(nèi)容更改）
showCategory =10, #只顯示前10
split="ONTOLOGY") + #以O(shè)NTOLOGY類型分開
facet_grid(ONTOLOGY~., scale='free') #以O(shè)NTOLOGY類型分開繪圖
dev.off()
此時在工作文件夾中得到了pdf格式的GO富集繪圖：
[圖片上傳失敗...(image-63c1b4-1679977340335)]
<figcaption style="color: rgb(153, 153, 153); font-size: 0.9em; line-height: 1.5; margin-top: calc(0.666667em); padding: 0px 1em; text-align: center;">barplot</figcaption>
3.2 dotplot
dotplot(eG,x = "GeneRatio", color = "p.adjust", size = "Count", #默認參數(shù)
showCategory =5,#只顯示前5
split="ONTOLOGY") + #以O(shè)NTOLOGY類型分開
facet_grid(ONTOLOGY~., scale='free') #以O(shè)NTOLOGY類型分屏繪圖
此時得到GO富集繪圖：
[圖片上傳失敗...(image-ff5536-1679977340335)]
<figcaption style="color: rgb(153, 153, 153); font-size: 0.9em; line-height: 1.5; margin-top: calc(0.666667em); padding: 0px 1em; text-align: center;">dotplot</figcaption>
4 ggplot2包繪圖
上面的圖有點單調(diào)妻导，我們可用ggplot2繪制更美觀一點的圖逛绵。
有些繪圖用的Enrichment Factor或者Fold Enrichment值為橫坐標。這個值需要自行計算倔韭。
Enrichment Factor = GeneRatio/BgRatio
GeneRatio：基因比术浪，分子是富集到此GO term上的基因數(shù)，而分母是所有得輸入基因數(shù)寿酌。
BgRatio：背景比胰苏，分子是此GO term得基因數(shù)，分母則是所有被GO注釋的基因數(shù)醇疼。
參考：生信交流平臺：GO和KEGG富集倍數(shù)（Fold Enrichment）如何計算
4.1 計算Enrichment Factor
本人用tidyverse包的separate函數(shù)將GeneRatio和BgRatio的分子分母先分開硕并，再進行計算法焰。
library(tidyverse)
eGo = read.table("eG.txt",header=TRUE,sep="\t",quote = "") #讀取第1部分enrichGO分析輸出的文件eG。
eGo <- separate(data=eGo, col=GeneRatio,into = c("GR1", "GR2"), sep = "/") #劈分GeneRatio為2列（GR1倔毙、GR2）
eGo <- separate(data=eGo, col=BgRatio, into = c("BR1", "BR2"), sep = "/") #劈分BgRatio為2列（BR1埃仪、BR2）
eGo <- mutate(eGo, enrichment_factor = (as.numeric(GR1)/as.numeric(GR2))/(as.numeric(BR1)/as.numeric(BR2))) #計算Enrichment Factor
此時eGo文件中已多了enrichment_factor列。
4.2 BP\MF\CC各取排名前10的term
只用排名前10的term畫圖陕赃，先把BP卵蛉、MF、CC三種類型各取top10么库，然后再組合在一起傻丝。
eGoBP <- eGo %>%
filter(ONTOLOGY=="BP") %>%
filter(row_number() >= 1,row_number() <= 10)
eGoCC <- eGo %>%
filter(ONTOLOGY=="CC") %>%
filter(row_number() >= 1,row_number() <= 10)
eGoMF <- eGo %>%
filter(ONTOLOGY=="MF") %>%
filter(row_number() >= 1,row_number() <= 10)
eGo10 <- rbind(eGoBP,eGoMF,eGoCC)
4.3 ggplot2畫圖
library(ggplot2)
p <- ggplot(eGo10,aes(enrichment_factor,Description)) +
geom_point(aes(size=Count,color=qvalue,shape=ONTOLOGY)) +
scale_color_gradient(low="red",high = "green") +
labs(color="pvalue",size="Count", shape="Ontology",
x="Enrichment Factor",y="GO term",title="GO enrichment") +
theme_bw()
p
[圖片上傳中...(image-a817a2-1679977340335-17)]
<figcaption style="color: rgb(153, 153, 153); font-size: 0.9em; line-height: 1.5; margin-top: calc(0.666667em); padding: 0px 1em; text-align: center;">ggplot2圖1</figcaption>
pvalue區(qū)分顏色，Count值區(qū)分大小诉儒，Ontology類型以不同形狀圖形表示葡缰。
4.4 ggplot2美化
由于p值都比較小，所以上圖顏色的區(qū)分度不夠允睹，可以用-log10(pvalue)值區(qū)分顏色运准。
p <- ggplot(eGo10,aes(enrichment_factor,Description)) +
geom_point(aes(size=Count,color=-1*log10(pvalue),shape=ONTOLOGY)) +
scale_color_gradient(low="green",high ="red") +
labs(color=expression(-log[10](p_value)),size="Count", shape="Ontology",
x="Enrichment Factor",y="GO term",title="GO enrichment") +
theme_bw()
p
[圖片上傳中...(image-5b4432-1679977340335-16)]
<figcaption style="color: rgb(153, 153, 153); font-size: 0.9em; line-height: 1.5; margin-top: calc(0.666667em); padding: 0px 1em; text-align: center;">ggplot2圖2</figcaption>
4.5 ggplot2分屏繪圖
以O(shè)NTOLOGY類型橫向分屏：
p + facet_wrap( ~ ONTOLOGY)
[圖片上傳中...(image-841e17-1679977340335-15)]
<figcaption style="color: rgb(153, 153, 153); font-size: 0.9em; line-height: 1.5; margin-top: calc(0.666667em); padding: 0px 1em; text-align: center;">ggplot2圖3</figcaption>
以O(shè)NTOLOGY類型縱向分屏：
p + facet_wrap( ~ ONTOLOGY,ncol= 1,scale='free')
[圖片上傳中...(image-9cc210-1679977340335-14)]
<figcaption style="color: rgb(153, 153, 153); font-size: 0.9em; line-height: 1.5; margin-top: calc(0.666667em); padding: 0px 1em; text-align: center;">ggplot2圖4</figcaption>
4.5 ggplot2調(diào)整排序
可以用fct_reorder(factor(x), y, .fun = median, .desc = FALSE)函數(shù)（將x按照y的順序排序）對繪圖排序。
參考：R語言學堂：forcats | fct_reorder2函數(shù)功能詳解及其在可視化中的應(yīng)用
這里將y軸Description按照x軸enrichment_factor的大小進行排序缭受。
aes(enrichment_factor, Description)要改成：
aes(enrichment_factor, fct_reorder(factor(Description), enrichment_factor))
p <- ggplot(eGo10,aes(enrichment_factor, fct_reorder(factor(Description), enrichment_factor))) +
geom_point(aes(size=Count,color=-1*log10(pvalue),shape=ONTOLOGY)) +
scale_color_gradient(low="green",high ="red") +
labs(color=expression(-log[10](p_value)),size="Count", shape="Ontology",
x="Enrichment Factor",y="GO term",title="GO enrichment") +
theme_bw()
p + facet_wrap( ~ ONTOLOGY,ncol= 1,scale='free')
[圖片上傳中...(image-ea471d-1679977340335-13)]
<figcaption style="color: rgb(153, 153, 153); font-size: 0.9em; line-height: 1.5; margin-top: calc(0.666667em); padding: 0px 1em; text-align: center;">ggplot2圖5</figcaption>
5 GOBubble圖
需要GOplot包
library(GOplot)
畫圖需要一個circle_dat(terms, genes)數(shù)據(jù)框
terms：包含 'category', 'ID', 'term', adjusted p-value ('adj_pval') 及'genes'的數(shù)據(jù)框胁澳。
genes：包含 'ID', 'logFC'的數(shù)據(jù)框。
5.1 準備circle_dat數(shù)據(jù)
這里我們使用eGo10數(shù)據(jù)（見4.2）米者，也就是前面BP韭畸、MF、CC三種類型各取top10的數(shù)據(jù)蔓搞。
GOterms = data.frame(category = eGo10$ONTOLOGY,
ID = eGo10$ID,
term = eGo10$Description,
genes = gsub("/", ",", eGo10$geneID),
adj_pval = eGo10$p.adjust)
[圖片上傳中...(image-672673-1679977340335-12)]
<figcaption style="color: rgb(153, 153, 153); font-size: 0.9em; line-height: 1.5; margin-top: calc(0.666667em); padding: 0px 1em; text-align: center;">GOterms部分內(nèi)容</figcaption>
genelist <- data.frame(ID = 數(shù)據(jù)$ID, logFC = 數(shù)據(jù)$logFC) #從已有“數(shù)據(jù)”中提取genelist胰丁，1列ID，1列l(wèi)ogFC喂分。
[圖片上傳中...(image-562051-1679977340335-11)]
<figcaption style="color: rgb(153, 153, 153); font-size: 0.9em; line-height: 1.5; margin-top: calc(0.666667em); padding: 0px 1em; text-align: center;">genelist格式</figcaption>
準備circle_dat數(shù)據(jù)
circ <- circle_dat(GOterms, genelist)
[圖片上傳中...(image-e37662-1679977340335-10)]
<figcaption style="color: rgb(153, 153, 153); font-size: 0.9em; line-height: 1.5; margin-top: calc(0.666667em); padding: 0px 1em; text-align: center;">circle_dat數(shù)據(jù)</figcaption>
5.2 繪制GOBubble圖
GOBubble(circ, labels = 5, # 標注的界值：-log(adjusted p-value) (默認5)
table.legend =T, #是否顯示右側(cè)legend锦庸，默認是
ID=T, # T標注term名稱，F(xiàn)標注ID
display='single') #是否分屏
[圖片上傳中...(image-29fd24-1679977340335-9)]
<figcaption style="color: rgb(153, 153, 153); font-size: 0.9em; line-height: 1.5; margin-top: calc(0.666667em); padding: 0px 1em; text-align: center;">GOBubble圖</figcaption>
GOBubble(circ, labels = 5, # 標注的界值：-log(adjusted p-value) (默認5)
table.legend =F, #不顯示右側(cè)legend
ID=F, # 標注term名稱
display='single') # 不分屏
[圖片上傳中...(image-89f8f4-1679977340335-8)]
<figcaption style="color: rgb(153, 153, 153); font-size: 0.9em; line-height: 1.5; margin-top: calc(0.666667em); padding: 0px 1em; text-align: center;">GOBubble圖2</figcaption>
<u style="border-bottom: 1px solid rgb(68, 68, 68); text-decoration: none;">z-score相當于上調(diào)基因數(shù)和下調(diào)基因數(shù)二者標準化的一個趨勢</u>
6 GOCircle圖
還是GOplot包蒲祈，也是用上面的circ數(shù)據(jù)甘萧，默認參數(shù)畫圖：
GOCircle(circ）
[圖片上傳中...(image-9d7a66-1679977340335-7)]
<figcaption style="color: rgb(153, 153, 153); font-size: 0.9em; line-height: 1.5; margin-top: calc(0.666667em); padding: 0px 1em; text-align: center;">GOCircle圖</figcaption>
參數(shù)設(shè)置：
GOCircle(circ,rad1=2, #內(nèi)環(huán)半徑
rad2=3, #外環(huán)半徑
label.size= 5, #標簽大小
label.fontface = 'bold', #標簽字體
nsub=10, #顯示的terms數(shù)，前10個梆掸。（畫圖前需要先把數(shù)據(jù)整理好扬卷，想要哪些term）
zsc.col = c('red', 'white', "green"), # z-score顏色
lfc.col = c('red', 'green')) # 基因up或down的顏色
[圖片上傳中...(image-520530-1679977340335-6)]
<figcaption style="color: rgb(153, 153, 153); font-size: 0.9em; line-height: 1.5; margin-top: calc(0.666667em); padding: 0px 1em; text-align: center;">GOCircle圖2</figcaption>
7 GOHeat熱圖
還是GOplot包，需要準備chord_dat(data, genes, process)數(shù)據(jù)：
le data-draft-node="block" data-draft-type="table" data-size="normal" data-row-style="normal">酸钦。
chord <- chord_dat(circ, #前面的circ數(shù)據(jù)
genelist[1:100,], #選擇需要的基因
GOterms$term[1:10]) #選擇要顯示的term
[圖片上傳中...(image-e73984-1679977340335-5)]
<figcaption style="color: rgb(153, 153, 153); font-size: 0.9em; line-height: 1.5; margin-top: calc(0.666667em); padding: 0px 1em; text-align: center;">chord_dat</figcaption>
GOHeat(chord, nlfc =1, #數(shù)據(jù)中有l(wèi)ogFC填1怪得，沒有填0，默認為0
fill.col = c('red', 'white', 'blue')) #自定義顏色
[圖片上傳中...(image-e73a6e-1679977340335-4)]
<figcaption style="color: rgb(153, 153, 153); font-size: 0.9em; line-height: 1.5; margin-top: calc(0.666667em); padding: 0px 1em; text-align: center;">GOHeat圖</figcaption>
8 GOChord弦圖
還是GOplot包，用chord_dat(data, genes, process)數(shù)據(jù)徒恋。
GOChord(chord)
[圖片上傳中...(image-337466-1679977340335-3)]
<figcaption style="color: rgb(153, 153, 153); font-size: 0.9em; line-height: 1.5; margin-top: calc(0.666667em); padding: 0px 1em; text-align: center;">GOChord弦圖</figcaption>
更改參數(shù)：
GOChord(chord, space = 0, #弦的間隔
gene.order = 'logFC', #基因排序方式
gene.space = 0.25, #基因名稱和圖的間隔
gene.size = 5, #基因名的大小
nlfc =1, #是否有l(wèi)ogFC
border.size= NULL, #彩虹邊框大小
lfc.col=c('red','black','cyan')) #自定義logFC 顏色
[圖片上傳中...(image-b81360-1679977340335-2)]
<figcaption style="color: rgb(153, 153, 153); font-size: 0.9em; line-height: 1.5; margin-top: calc(0.666667em); padding: 0px 1em; text-align: center;">GOChord弦圖2</figcaption>
9 樹狀圖
eGoBP <- enrichGO(gene = g_list,
OrgDb = org.Hs.eg.db,
pvalueCutoff =0.01,
qvalueCutoff = 0.01,
ont="BP", #BP\MF\CC
readable =T)
plotGOgraph(eGoBP)
[圖片上傳中...(image-ce8de8-1679977340334-1)]
<figcaption style="color: rgb(153, 153, 153); font-size: 0.9em; line-height: 1.5; margin-top: calc(0.666667em); padding: 0px 1em; text-align: center;">plotGOgraph</figcaption>
