好久没更新了
今天就更新一波，今天的主题是如何使用ggplot2个性化绘制GO富集分析柱状图。
进入正题,为什么选择利用ggplot做图，而没用clusterProfiler自带的barplot或dotplot作图，利用ggplot最好的一个是更容易个性化选择自己想要展示的GO条目，这往往为我们讲好自己的生物学故事有很大的帮助。由于每个人需要讲述的生物学故事不同，我这里就只是展示最常规的做法，选取FDR值最小的Top10条目进行展示。
1.使用clusterProfiler包进行GO富集分析
##上下调基因的GO功能分析
#ID转换
library(clusterProfiler)
library(org.Hs.eg.db)
up.genes = data$gene[which(data$Group == "up-regulated")]
up.genes.df <- bitr(up.genes, fromType = "SYMBOL", #fromType是指你的数据ID类型是属于哪一类的
toType = c("ENSEMBL", "ENTREZID"), #toType是指你要转换成哪种ID类型，可以写多种，也可以只写一种
OrgDb = org.Hs.eg.db)#Orgdb是指对应的注释包是哪个
head(up.genes.df)
down.genes = data$gene[which(data$Group == "down-regulated")]
down.genes.df <- bitr(down.genes, fromType = "SYMBOL", #fromType是指你的数据ID类型是属于哪一类的
toType = c("ENSEMBL", "ENTREZID"), #toType是指你要转换成哪种ID类型，可以写多种，也可以只写一种
OrgDb = org.Hs.eg.db)#Orgdb是指对应的注释包是哪个
head(down.genes.df)
#上调基因GO富集
up.gene=up.genes.df$ENTREZID
#GO富集分析
up.kk <- enrichGO(gene = up.gene,
OrgDb = org.Hs.eg.db,
pvalueCutoff =0.05,
qvalueCutoff = 0.05,
ont="all",
readable =T)
#下调基因GO富集
down.gene=down.genes.df$ENTREZID
#GO富集分析
down.kk <- enrichGO(gene = down.gene,
OrgDb = org.Hs.eg.db,
pvalueCutoff =0.05,
qvalueCutoff = 0.05,
ont="all",
readable =T)
2.ggplot2绘制GO条目图（我只关注BP）
2.1 按照qvalue升序排序，分别选出前10个BP,CC,MF的条目（只展示Top10 BP条目），由于enrichGO函数生成的数据框默认是按照qvalue升序排序，所以这里我们只用选取前十个就行了
先画up-regulated genes GO terms
#go_MF<-go[go$ONTOLOGY=="MF",][1:10,]
#go_CC<-go[go$ONTOLOGY=="CC",][1:10,]
go_BP<-up.kk[up.kk$ONTOLOGY=="BP",][1:10,]
go_enrich_df<-data.frame(ID=go_BP$ID,
Description=go_BP$Description,
GeneNumber=go_BP$Count,
type=factor(rep("biological process", 10),levels="biological process"))
## numbers as data on x axis
go_enrich_df$number <- factor(rev(1:nrow(go_enrich_df)))
## add GO terms description names
labels <- go_enrich_df$Description
names(labels) = rev(1:nrow(go_enrich_df))
## colors for bar // green, blue, orange
CPCOLS <- "#eb6bac"
## optional other colors
#CPCOLS <- c("#8DA1CB", "#FD8D62", "#66C3A5")
library(ggplot2)
p1 <- ggplot(data=go_enrich_df, aes(x=number, y=GeneNumber, fill=type)) +
geom_bar(stat="identity", width=0.8) + coord_flip() +
scale_fill_manual(values = CPCOLS) + theme_test() +
scale_x_discrete(labels=labels) +
xlab("GO term") +
theme(axis.text=element_text(face = "bold", color="gray50")) +
labs(title = "The Most Enriched GO Terms for Up.genes")#coord_flip(...)横向转换坐标：把x轴和y轴互换，没有特殊参数
ggsave("01.up.gene_GO_BP.pdf", width = 10, height = 8)
01.up.gene_GO_BP.png
接着画down-regulated genes GO terms
go_BP<-down.kk[down.kk$ONTOLOGY=="BP",][1:10,]
go_enrich_df<-data.frame(ID=go_BP$ID,
Description=go_BP$Description,
GeneNumber=go_BP$Count,
type=factor(rep("biological process", 10),levels="biological process"))
## numbers as data on x axis
go_enrich_df$number <- factor(rev(1:nrow(go_enrich_df)))
## add GO terms description names
labels <- go_enrich_df$Description
names(labels) = rev(1:nrow(go_enrich_df))
## colors for bar // green, blue, orange
CPCOLS <- "#ffa8a0"
## optional other colors
#CPCOLS <- c("#8DA1CB", "#FD8D62", "#66C3A5")
library(ggplot2)
p2 <- ggplot(data=go_enrich_df, aes(x=number, y=GeneNumber, fill=type)) +
geom_bar(stat="identity", width=0.8) + coord_flip() +
scale_fill_manual(values = CPCOLS) + theme_test() +
scale_x_discrete(labels=labels) +
xlab("GO term") +
theme(axis.text=element_text(face = "bold", color="gray50")) +
labs(title = "The Most Enriched GO Terms for Down.genes")#coord_flip(...)横向转换坐标：把x轴和y轴互换，没有特殊参数
ggsave("02.down.gene_GO_BP.pdf", width = 10, height = 8)
02.down.gene_GO_BP.png
最后合并出图
library(ggpubr)#利用ggarrange将图像排列在一起
ggarrange(p1/p2,widths = c(1,1),heights = c(1,1))
ggsave("03.GO_BP.pdf", width = 12, height = 8)
03.GO_BP.png
中午花了睡觉的时间写完啦，下午去放松玩去了，哈哈哈哈
