GO富集
GO总共有三个ontology(本体),分别描述基因的分子功能(molecular function)、细胞组分(cellular component)、参与的生物过程(biological process)。GO的基本单位是term(词条、节点),每个term都对应一个属性。
准备你感兴趣的基因列表文件示例如下
开始分析部分
rm(list=ls())
#可以设置一下工作路径
setwd("D:/RStudio/R-work/GO条形图添加基因数量标签test")
getwd()#设置完检查一下是否在这个路径中
#加载需要的包
library(clusterProfiler)#富集分析包
library(org.Hs.eg.db)#指定物种,我用的人的。
library(ggplot2)#我加载这个报的目的是给柱形图加基因数量标签
#基因list文件
geneinfo <- read.table("gene2.txt", header = T, sep = '\t')
geneinfo$ID <- as.character(geneinfo$ID) #character格式
# 将SYMBOL格式转为ENSEMBL和ENTERZID格式 ,也可以是其中一种
genelist <- bitr(geneinfo$ID, fromType="SYMBOL", toType=c("ENSEMBL", "ENTREZID"), OrgDb="org.Hs.eg.db")#在这里遇到报错,如图“转换格式报错”这里是因为我选择的测试基因在物种中一个都没有比对到(仅代表个人想法)
head(genelist,2)#查看一下转换结果
write.table(genelist,"genelist.csv", sep=",") #可以将genelist数据导出
# enrichGO 富集分析
go_enrich <- enrichGO(gene = genelist$ENTREZID,
OrgDb = org.Hs.eg.db, #没有organism="human",改为OrgDb=org.Hs.eg.db
#keytype = 'ENSEMBL',
ont = "ALL", #也可以是 CC BP MF中的一种
pAdjustMethod = "BH", #矫正方式 holm”, “hochberg”, “hommel”, “bonferroni”, “BH”, “BY”, “fdr”, “none”中的一种
pvalueCutoff = 1, #P值会过滤掉很多,可以全部输出
qvalueCutoff = 1,
readable = TRUE)
head(go_enrich,2)
# 输出结果
write.csv(summary(go_enrich),"goenrich.csv",row.names =FALSE)
# 绘图
##可视化--点图
dotplot(go_enrich,title="EnrichmentGOdot",font.size = 15)#点图,按富集的数从大到小的
##########可视化--条形图
#加数字标签的
barplot(go_enrich, showCategory=20,font.size=10,title="EnrichmentGO",label_format=100)+#条状图,按p从小到大排,绘制前20个Term
geom_text(aes(label=Count),family="sans",position = position_stack(1.05),size=3)
#设置了左侧条目长度label_format
barplot(go_enrich, font.size=10,title="EnrichmentGO",label_format=100)
#将富集分析三个部分分别展示10个条目
barplot(go_enrich, split="ONTOLOGY",font.size=10,showCategory = 10,label_format=100) + facet_grid(ONTOLOGY~., scale="free")
########如果想要将图片输出(这个尺寸不合适,需要再次调整)
pdf("barplot2.pdf",width=50,height=100)
barplot(go_enrich, split="ONTOLOGY",showCategory = 10,label_format=100)
dev.off()
##可视化--
plotGOgraph(go_enrich)#这个应该是BP\MF\CC其中一个,才可以出这个图