为啥搞这个
之前写过一个用ggplot2画TBtools富集分析结果的代码，用ggplot2做富集分析气泡图, 发现大家确实喜欢花里胡哨的东西.... 其实我感觉柱状图就已经挺能说明问题的，当然一张花里胡哨的图小概率可以成为论文的加分项，这里就仿照Y叔 clusterProfiler里面两个富集合并出图的结果用ggplot2实现了下。
代码如下
最近比较懒，没有打getopt，打开Rstudio跑下吧，或者保存成.R文件，到时候source一下也可以
godotplot = function(left,right,lname,rname,name){
## top 15 terms
x = left[order(left$P_value)[1:15],]
y = right[order(right$P_value)[1:15],]
## Gene ratio calculate
x$GeneRatio = x$HitsGenesCountsInSelectedSet/x$AllGenesCountsInSelectedSet
y$GeneRatio = y$HitsGenesCountsInSelectedSet/y$AllGenesCountsInSelectedSet
## ordered the GO name for plot
x = x[order(x$GeneRatio),]
y = y[order(y$GeneRatio,decreasing = T),]
## data cleanning
leftdata = data.frame(GO_Name = factor(x$GO_Name,levels = x$GO_Name),
P_value = -log10(x$P_value),
GeneRatio = x$GeneRatio,
Cultivar = rname)
rightdata = data.frame(GO_Name = factor(y$GO_Name,levels = y$GO_Name),
P_value = -log10(y$P_value),
GeneRatio = y$GeneRatio,
Cultivar = lname)
## merge
z = rbind(leftdata,rightdata)
## plot
p = ggplot(data = z, mapping = aes(x = Cultivar,y = GO_Name))+
geom_point(aes(size = GeneRatio,color = P_value))+
scale_color_gradient(low = "blue", high = "red", guide=guide_colorbar(reverse=TRUE))+
scale_y_discrete(labels=function(y) stringr::str_wrap(y,width=50)) +
theme_bw() +
theme(axis.text.x = element_text(colour = "black",
size = 16, vjust =1 ),
axis.text.y = element_text(colour = "black",
size = 16, hjust =1 ),
axis.title = element_text(margin=margin(10, 5, 0, 0),
color = "black",size = 18),
axis.title.y = element_text(angle=90))+
labs(color=expression(-log[10](Qvalue)),
size="Gene Ratio",
x = "",
y = "")
return(p)
# ggsave(p,filename = paste("~/yourpath",name,".pdf",sep = ""),width = 9,height = 10,dpi = 300)
}
说明
- 我设置了按照
P_value排序的前15个，如果要改规则，在function里面改一下 - Y轴排序是按照
GeneRatio排的，想改规则也对应着修改一下 - 如果要直接保存的话把return(p)注释掉，去掉ggsave前面的#
- 其实你读懂代码的话可以合并无限个.....
用法
这个就比较简单了，
先用TBtools撸两个富集分析的结果出来，然后保存好。
## step1 import TBtools enrichment result
left = read.table("left.xls", header = T, sep ="\t")
right = read.table("right.xls", header = T, sep = "\t")
## if return(p) was comment but ggsave(...) not
godotplot(left = left,right = right,lname = "left",rname = "right",name = "testGOdotplot")
## if ggsave(...) was comment but return(p) not
p = godotplot(left = left,right = right,lname = "left",rname = "right",name = "testGOdotplot")
p
ggsave(p,"filepath/filename", width = 9, height = 10,dpi = 300)
完了 最后看下效果图吧：goterm被截掉了哈..
testGOdotplot.png
