已知文件如下
原始文件
- DEG.list #差异基因list
- rice.map #物种注释到GO数据库的信息 一般只用到上面两个文件就可以
- rice.map # 1. blast2go 软件做注释，再使用 mapping2maplist.pl 脚本进行格式转换 # 2. 可以直接从数据库下载，再使用 plaza2maolist.pl 脚本进行格式转化
原始数据上传到百度网盘：
链接：https://pan.baidu.com/s/1QVd-s9PWY9f0LpMCQcsyFw
提取码：51og
下面放出完整操作R代码，只需要修改导入的文件就可直接运行得到结果
##############################################
##### 2020/3/6 w
##### topGO
##############################################
# 安装topGO软件包
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("topGO", version = "3.8")
BiocManager::install("Rgraphviz", version = "3.8")
# 设置工作目录，后面读取文件什么的就可以直接读取不需要那么长的路径
setwd('D:/test_data')
# 加载包
rm(list=ls())
library(topGO)
library(Rgraphviz)
# 设置输入文件，后面直接在这个地方修改文件名称就可以直接运行了
input="DEG.list" #差异基因名称的列表
mapfile="rice.map" #所有基因GO map结果，也可以是 list 文件等
# 开始分析
gene_id = readMappings(file = mapfile) #如果是读取其他文件格式，后面参数还需要修改
gene_names = names(gene_id)
my_genes = read.table(input)[,1]
gene_list = rep(1,length(gene_id))
names(gene_list) = names(gene_id)
gene_list[match(my_genes,names(gene_list))] = 0
top_diff_genes = function(allScore){return(allScore<0.01)}
# BP 富集分析
#new() 创建一个 topGO 的对象，然后对这个对象做检验
sample_gOdata = new("topGOdata",
nodeSize = 6,
ontology="BP",
allGenes = gene_list,
annot = annFUN.gene2GO,
gene2GO = gene_id,
geneSel=top_diff_genes)
# 做检验，使用的是elim 的算法，使用 ks 的统计量。可以理解为 p 值
result_KS.elim = runTest(sample_gOdata,
algorithm = "elim",
statistic = "ks")
#提取基因 table
allres = GenTable(sample_gOdata,
KS = result_KS.elim,
ranksOf = "classic",
topNodes = attributes(result_KS.elim)$geneData[4])
#生成文件，后面画图都可以用这个表
write.table(allres,
file = paste(input,".BP.xls",sep=""),
sep="\t", quote=FALSE, col.names=TRUE, row.names=FALSE)
# 输出矢量图
pdf(paste(input,".BP.pdf",sep=""))
showSigOfNodes(sample_gOdata,
score(result_KS.elim),
firstSigNodes = 10,
useInfo = "all") #设置节点数量，10个或者20个更多都可以
dev.off()
# 输出像素图
png(paste(input,".BP.png",sep=""))
showSigOfNodes(sample_gOdata, score(result_KS.elim), firstSigNodes = 10, useInfo = "all")
dev.off()
# MF 富集分析（同理）
sample_gOdata = new("topGOdata",
nodeSize = 6,
ontology="MF",
allGenes = gene_list,
annot = annFUN.gene2GO,
gene2GO = gene_id,
geneSel=top_diff_genes)
# 做检验，使用的是elim 的算法，使用 ks 的统计量。可以理解为 p 值
result_KS.elim = runTest(sample_gOdata,
algorithm = "elim",
statistic = "ks")
#提取基因 table
allres = GenTable(sample_gOdata,
KS = result_KS.elim,
ranksOf = "classic",
topNodes = attributes(result_KS.elim)$geneData[4])
#生成文件，后面画图都可以用这个表
write.table(allres,
file = paste(input,".MF.xls",sep=""),
sep="\t", quote=FALSE, col.names=TRUE, row.names=FALSE)
# 输出矢量图
pdf(paste(input,".MF.pdf",sep=""))
showSigOfNodes(sample_gOdata,
score(result_KS.elim),
firstSigNodes = 10,
useInfo = "all") #设置节点数量，10个或者20个更多都可以
dev.off()
# 输出像素图
png(paste(input,".MF.png",sep=""))
showSigOfNodes(sample_gOdata, score(result_KS.elim), firstSigNodes = 10, useInfo = "all")
dev.off()
#CC节点的富集分析
sample_gOdata = new("topGOdata",
nodeSize = 6,
ontology="CC",
allGenes = gene_list,
annot = annFUN.gene2GO,
gene2GO = gene_id,
geneSel=top_diff_genes)
# 做检验，使用的是elim 的算法，使用 ks 的统计量。可以理解为 p 值
result_KS.elim = runTest(sample_gOdata,
algorithm = "elim",
statistic = "ks")
#提取基因 table
allres = GenTable(sample_gOdata,
KS = result_KS.elim,
ranksOf = "classic",
topNodes = attributes(result_KS.elim)$geneData[4])
#生成文件，后面画图都可以用这个表
write.table(allres,
file = paste(input,".CC.xls",sep=""),
sep="\t", quote=FALSE, col.names=TRUE, row.names=FALSE)
# 输出矢量图
pdf(paste(input,".CC.pdf",sep=""))
showSigOfNodes(sample_gOdata,
score(result_KS.elim),
firstSigNodes = 10,
useInfo = "all") #设置节点数量，10个或者20个更多都可以
dev.off()
# 输出像素图
png(paste(input,".CC.png",sep=""))
showSigOfNodes(sample_gOdata, score(result_KS.elim), firstSigNodes = 10, useInfo = "all")
dev.off()
输出table文件
- GO.ID # 富集到的GO.ID
- Term # 功能的描述
- Annotated # 该节点注释到的基因数目
- Significant # 差异基因有多少个注释到该节点
- Expected # 期望有多少个
- KS # 相当于 p 值，越小越显著
可以用输出文件做出更好的可视化结果
- 本文参考文件
腾讯课堂视频 PlantTech（普朗泰科）学院 —— R 应用- topGO 富集分析
topGO文档
链接：https://pan.baidu.com/s/1OsqGVDcD1XZtpj5BTgncZg
提取码：08r8
