栏目分类:
子分类:
返回
名师互学网用户登录
快速导航关闭
当前搜索
当前分类
子分类
实用工具
热门搜索
名师互学网 > IT > 软件开发 > 后端开发 > Go语言

R语言|GO富集分析

Go语言 更新时间: 发布时间: IT归档 最新发布 模块sitemap 名妆网 法律咨询 聚返吧 英语巴士网 伯小乐 网商动力

R语言|GO富集分析

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

一、数据链接
链接: https://pan.baidu.com/s/1PPUW5YyJHjwxvkjmMi0Vtw?pwd=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="t",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="t",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) 

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

转载请注明:文章转载自 www.mshxw.com
本文地址:https://www.mshxw.com/it/902410.html
我们一直用心在做
关于我们 文章归档 网站地图 联系我们

版权所有 (c)2021-2022 MSHXW.COM

ICP备案号:晋ICP备2021003244-6号