在进行富集分析绘图的时候,有的时候一次需要绘制多张GO和KEGG图片,可以通过简单的脚本实现批量绘制。
1. 所需R包
(1)ggplot2
(2)clusterProfiler
(3)AnnotationHub
2.分析所用数据
大豆根系不同时期的WGCNA分析,对获得的共表达数据模块进行富集分析(数据已经完成了ID转换和GO富集)。
# 1. 安装所需要的R包,如果已经安装不需要再次安装。
install.packages("ggplot2")
library("BiocManager")
BiocManager::install(c("clusterProfiler","AnnotationHub"))
# 2. 安装所需的注释库
# 构建hub
hub <- AnnotationHub::AnnotationHub()
query(hub,"Glycine max") # 查找其中的大豆数据
Glycine.db = hub[["AH85411"]]
# 编写R脚本
#!/usr/bin/R
suppressPackageStartupMessages(library(clusterProfiler))
suppressPackageStartupMessages(library(AnnotationHub))
suppressPackageStartupMessages(library(ggplot2))
# 输入和输出文件
args <- commandArgs(trailingOnly=TRUE)
inFile <- args[1]
outFile1 <- paste0(strsplit(inFile,split="\\.")[[1]][3],".barplot.",".png")
outFile2 <- paste0(strsplit(inFile,split="\\.")[[1]][3],".Dotplot.",".png")
# 有的富集分析的描述信息太长,适当的缩短
# 例如:zinc finger AN1 and C2H2 domain-containing stress-associated protein 16
cutStr <- function(x){
if(nchar(x)>=30){
x = strsplit(x,split=" ")[[1]][1:min(length(strsplit(x," ")),6)] # 描述信息,字符串长度大于30,截取前面的6个字符串。
x = paste0(x,"...")
}else{
return(x)
}
}
# 读入文件(文件使用clusterProfiler已经完成了Gene ID转换)
df <- read.table(inFile,header=TRUE,stringsAsFactors=FALSE,sep="\t")
# ev <- function(x){eval(parse(text=x))}
df$generation <- round(sapply(df$GeneRatio,function(x){eval(parse(text=x))}),3) ## 将GeneRation转换为小数
df$new_Description <- sapply(df$Description,cutStr)
# 开始绘图
## 绘制条型图
p1 <- ggplot(data=df)+ geom_bar(aes(x=new_Description,y=-log10(pvalue), fill=ONTOLOGY), stat='identity')+
coord_flip() + scale_x_discrete(limits=df$new_Description)+labs(title=strsplit(inFile,split="\\.")[[1]][1])+
theme(plot.title = element_text(hjust=0.5))
ggsave(file=outFile1,p1,width = 8, height=7)
## 绘制点图
p2 <- ggplot(df,aes(x=generation,y=new_Description,color=-1*log10(qvalue)))+
geom_point(aes(size=Count))+
labs(title=strsplit(inFile,split="\\.")[[1]][1])+
scale_color_gradient(low = 'blue', high = 'red')+
theme(plot.title = element_text(hjust=0.5))
ggsave(file=outFile2,p2,width = 8, height=7)
3. 在shell环境中调用脚本
# 将上面的脚本保存为,Get_GO_Profiler.R
ls *.GO.txt | while read id;do Rscript barplotGO.R;done
## 注释文件的数目如果很多的化,也可以使用shell并发。