R语言clusterProfiler包GO富集分析（enrichplot包、GOplot包和ggplot2绘图）

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

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>

最后编辑于：2023.03.30 18:23:31

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 203,547评论 6赞 477
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 85,399评论 2赞 381
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 150,428评论 0赞 337
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 54,599评论 1赞 274
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 63,612评论 5赞 365
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 48,577评论 1赞 281
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 37,941评论 3赞 395
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,603评论 0赞 258
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 40,852评论 1赞 297
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,605评论 2赞 321
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,693评论 1赞 329
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,375评论 4赞 318
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 38,955评论 3赞 307
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 29,936评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 31,172评论 1赞 259
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 43,970评论 2赞 349
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 42,414评论 2赞 342

R语言clusterProfiler包GO富集分析（enrichplot包、GOplot包和ggplot2绘图）

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 树状图

推荐阅读更多精彩内容