https://zhuanlan.zhihu.com/p/597338272
https://zhuanlan.zhihu.com/p/569459288
https://www.jianshu.com/p/25386890a751
https://www.jianshu.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函数富集分析
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调整排序
5 GOBubble图
5.1 准备circle_dat数据
5.2 绘制GOBubble图
6 GOCircle图
7 GOHeat热图
8 GOChord弦图
9 树状图
通过某些方法（如差异基因、相关基因等）获得一批基因后，可以进行功能富集分析（包括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函数富集分析
eG <- enrichGO(gene = g_list, #需要分析的基因的EntrezID
OrgDb = org.Hs.eg.db, #人基因数据库
pvalueCutoff =0.01, #设置pvalue界值
qvalueCutoff = 0.01, #设置qvalue界值(FDR校正后的p值）
ont="all", #选择功能富集的类型，可选BP、MF、CC，这里选择all。
readable =T)
富集分析的类型可选BP（biological process）、MF（molecular function）、CC（Cellular Component）。
输出结果为txt文件
write.table(eG,file="eG.txt", sep="\t", quote=F, row.names = F)
结果如下：
[图片上传失败...(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", #默认参数（x和color可以根据eG里面的内容更改）
showCategory =10, #只显示前10
split="ONTOLOGY") + #以ONTOLOGY类型分开
facet_grid(ONTOLOGY~., scale='free') #以ONTOLOGY类型分开绘图
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", #默认参数
showCategory =5,#只显示前5
split="ONTOLOGY") + #以ONTOLOGY类型分开
facet_grid(ONTOLOGY~., scale='free') #以ONTOLOGY类型分屏绘图
此时得到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包绘图
上面的图有点单调，我们可用ggplot2绘制更美观一点的图。
有些绘图用的Enrichment Factor或者Fold Enrichment值为横坐标。这个值需要自行计算。
Enrichment Factor = GeneRatio/BgRatio
GeneRatio：基因比，分子是富集到此GO term上的基因数，而分母是所有得输入基因数。
BgRatio：背景比，分子是此GO term得基因数，分母则是所有被GO注释的基因数。
参考：生信交流平台：GO和KEGG富集倍数（Fold Enrichment）如何计算
4.1 计算Enrichment Factor
本人用tidyverse包的separate函数将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区分颜色，Count值区分大小，Ontology类型以不同形状图形表示。
4.4 ggplot2美化
由于p值都比较小，所以上图颜色的区分度不够，可以用-log10(pvalue)值区分颜色。
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分屏绘图
以ONTOLOGY类型横向分屏：
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>
以ONTOLOGY类型纵向分屏：
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调整排序
可以用fct_reorder(factor(x), y, .fun = median, .desc = FALSE)函数（将x按照y的顺序排序）对绘图排序。
参考：R语言学堂：forcats | fct_reorder2函数功能详解及其在可视化中的应用
这里将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)数据框
terms：包含 'category', 'ID', 'term', adjusted p-value ('adj_pval') 及'genes'的数据框。
genes：包含 'ID', 'logFC'的数据框。
5.1 准备circle_dat数据
这里我们使用eGo10数据（见4.2），也就是前面BP、MF、CC三种类型各取top10的数据。
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部分内容</figcaption>
genelist <- data.frame(ID = 数据$ID, logFC = 数据$logFC) #从已有“数据”中提取genelist，1列ID，1列logFC。
[图片上传中...(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数据
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数据</figcaption>
5.2 绘制GOBubble图
GOBubble(circ, labels = 5, # 标注的界值：-log(adjusted p-value) (默认5)
table.legend =T, #是否显示右侧legend，默认是
ID=T, # T标注term名称，F标注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, #不显示右侧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相当于上调基因数和下调基因数二者标准化的一个趋势</u>
6 GOCircle图
还是GOplot包，也是用上面的circ数据，默认参数画图：
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>
参数设置：
GOCircle(circ,rad1=2, #内环半径
rad2=3, #外环半径
label.size= 5, #标签大小
label.fontface = 'bold', #标签字体
nsub=10, #显示的terms数，前10个。（画图前需要先把数据整理好，想要哪些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)数据：
le data-draft-node="block" data-draft-type="table" data-size="normal" data-row-style="normal">。
chord <- chord_dat(circ, #前面的circ数据
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, #数据中有logFC填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)数据。
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>
更改参数：
GOChord(chord, space = 0, #弦的间隔
gene.order = 'logFC', #基因排序方式
gene.space = 0.25, #基因名称和图的间隔
gene.size = 5, #基因名的大小
nlfc =1, #是否有logFC
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>
