GO富集分析 —— 差异圈图、聚类圈图

一、数据链接
链接: https://pan.baidu.com/s/1PPUW5YyJHjwxvkjmMi0Vtwpwd=c9s8
提取码: c9s8

二、安装包介绍
clusterProfiler包：功能也比较强大，主要是做GO和KEGG的功能富集及其可视化。
org.Hs.eg.db包：转换NCBI、ensemble等数据库中基因ID，symbol等之间的转换。
enrichplot包：实现多种可视化方法来解释富集结果。
GOplot：功能富集绘图。

三、绘图
1、安装加载包

#安装包
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("org.Hs.eg.db")
BiocManager::install("DOSE")
BiocManager::install("clusterProfiler")
BiocManager::install("enrichplot")
install.packages("colorspace")
install.packages("stringi")
install.packages("ggplot2")
install.packages("digest")
install.packages("GOplot")
#加载包
library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)
library(ggplot2)
library(stringi)	# 处理表格数据的包
library(GOplot)

2、设置工作路径

setwd("D:\demo\GOcircos")

3、数据整理

inputFile="input.txt" 
# 读取文件      
rt=read.table(inputFile,sep="	",header=T,check.names=F) 

genes=as.vector(rt[,1])	# 选取rt第一列gene保存到genes
entrezIDs=mget(genes, org.Hs.egSYMBOL2EG, ifnotfound=NA)	# 找出基因对应的ID
entrezIDs=as.character(entrezIDs)	# 获取数据
rt=cbind(rt,entrezID=entrezIDs)		# 添加一列entrezID
rt=rt[is.na(rt[,"entrezID"])==F,]	# 删除没有基因的ID
gene=rt$entrezID

GO富集分析

GO=enrichGO(gene = gene,
            OrgDb = org.Hs.eg.db, # 参考基因组
            pvalueCutoff =1,	# P值阈值
            qvalueCutoff = 1,	# qvalue是P值的校正值
            ont="all",	# 主要的分为三种，三个层面来阐述基因功能，生物学过程（BP），细胞组分（CC），分子功能（MF）
            readable =T)	# 是否将基因ID转换为基因名
# 强制转换为数据框
GO=as.data.frame(GO)

# 筛选显著富集的数据
GO<-GO[(GO$pvalue<0.05 & GO$p.adjust<0.05),]
# 保存数据
write.table(GO,file="GO1.txt",sep="	",quote=F,row.names = F)
# 构建数据框矩阵
go=data.frame(Category = GO$ONTOLOGY,ID = GO$ID,Term = GO$Description, Genes = gsub("/", ", ", GO$geneID), adj_pval = GO$p.adjust)
# 构建数据框矩阵
genelist=data.frame(ID = rt$gene, logFC = rt$logFC)
row.names(genelist)=genelist[,1]
circ <- circle_dat(go, genelist)
termNum = 5	#限定GO数目
termNum=ifelse(nrow(go)<termNum,nrow(go),termNum)
geneNum = nrow(genelist)	#限定基因数目

4、绘图、保存图片

# 差异圈图
chord <- chord_dat(circ, genelist[1:geneNum,], go$Term[1:termNum])
pdf(file="GOcircos.pdf",width = 10,height = 10.2)
GOChord(chord, 
        space = 0.001,           # 基因之间的间距
        gene.order = 'logFC',    # 排序基因
        gene.space = 0.25,       # 基因离圆圈距离
        gene.size = 3,           # 基因字体大小
        border.size = 0.1,		 # 线的大小
        process.label = 7)       # GO名称大小
dev.off()

#聚类圈图
pdf(file="GOcluster.pdf",width = 12,height = 9)
GOCluster(circ, as.character(go[1:termNum,3]))
dev.off()

END