R中使用topGO进行富集分析

已知文件如下

• 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()

• GO.ID # 富集到的GO.ID
• Term # 功能的描述
• Annotated # 该节点注释到的基因数目
• Significant # 差异基因有多少个注释到该节点
• Expected # 期望有多少个
• KS # 相当于 p 值，越小越显著

可以用输出文件做出更好的可视化结果

• 本文参考文件
  腾讯课堂视频 PlantTech（普朗泰科）学院 —— R 应用- topGO 富集分析
  topGO文档
  链接：https://pan.baidu.com/s/1OsqGVDcD1XZtpj5BTgncZg
  提取码：08r8