1 概述
GO (Gene Ontology) 富集分析是最常见的生信分析需求之一，与之类似的还有KEGG通路富集分析、基因家族分析等等，这里我们直奔主题，给大家分享一些常用的GO富集分析工具以及对应的使用方法。
2 工具推荐
2.1 clusterprofiler的安装和使用
clusterprofiler是一个发布在Bioconductor上用于组学富集分析的R包工具，目前已经更新到了第4个版本，有关clusterprofiler的详细使用方法，请参考官方的说明文档，这里只是为了帮助大家能够快速上手完成一个基因富集分析的任务。
安装
# 如果下载速度太慢，可以考虑使用镜像下载
# options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor")
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("clusterProfiler")
BiocManager::install("AnnotationForge")
数据准备
Gene2GOs.map - 目标物种的背景注释文件，总共两列，列与列中间用制表符（即Tab键）分隔，每一个GO ID (GO:xxxxxxx) 会对应一个基因 ID，格式如下：
gene_ID1 GO_ID1
gene_ID2 GO_ID2
gene_ID3 GO_ID3
...
GeneList.csv - 用于分析的目标基因文件，可选择只保留一列用于富集分析的基因ID，如果需要其余辅助信息，如fold change、pvalue等，最好将基因ID放在第一列，格式如下：
gene_ID1
gene_ID2
gene_ID3
...
简单使用方法
先制作目标物种的OrgDb数据对象
library(AnnotationForge)
# 清除变量信息
rm(list=ls())
# 重新设置你的工作目录，你可以把它设置为数据文件所在的目录
setwd(".")
# 这里只是演示如何制作一个最简单的OrgDb，想要制作更详细的OrgDb请参考官方教程
# 读取数据
geneID2GO <- read.table(file = 'Gene2GOs.map', header = FALSE)
colnames(geneID2GO) <- c("GID", "GO")
# 数据处理
# 基因信息，3列："GID","SYMBOL","GENENAME"
fSym <- transform(geneID2GO, SYMBOL = GID, GENENAME = GID)[, c("GID","SYMBOL","GENENAME")]
fSym <- fSym[fSym[,2]!="-",]
fSym <- fSym[fSym[,3]!="-",]
fSym <- fSym[!duplicated(fSym),]
dim(fSym)
# 基因与GO的映射关系，3列："GID", "GO"，"EVIDENCE"
# gene id, go id (like GO:xxxxxxx) and evidence codes (like ISO)
fGO <- transform(geneID2GO, EVIDENCE=rep("", length(GID)))
fGO <- fGO[!duplicated(fGO),]
dim(fGO)
# 制作OrgDb
# 这里以油菜为例，邮箱和作者信息可以替换为实际的名字
makeOrgPackage(gene_info=fSym, go=fGO,
version="0.2",
maintainer="lxd <xxxxxxxxx@163.com>",
author="lxd <xxxxxxxxx@163.com>",
outputDir = ".",
tax_id="3708",
genus="Brassica",
species="napus",
goTable="go")
# 把制作好的OrgDb安装到本地，用于后续的加载
install.packages("./org.Bnapus.eg.db", type="source", repos=NULL)
library(clusterProfiler)
library(org.Bnapus.eg.db)
# 清除变量信息
rm(list=ls())
# 重新设置你的工作目录，你可以把它设置为数据文件所在的目录
setwd(".")
# step1: 读取目标基因文件
genes <- read.csv("GeneList.csv", header = FALSE)
geneList <- genes[, 1]
# step2: 执行富集检验
ego <- enrichGO(gene = geneList,
OrgDb = org.Bnapus.eg.db,
ont = "BP",
keyType = 'GID',
pAdjustMethod = "BH",
pvalueCutoff = 1,
qvalueCutoff = 1,
minGSSize = 1)
head(ego)
# step3: 结果分析和可视化
## 保存结果
write.csv(ego, paste("output", "clusterProfiler.csv", sep = "_"), row.names = FALSE)
## Bar plot
p1 <- barplot(ego, showCategory=10)
ggsave("GO_barplot.png", width=7, height=7)
ggsave("GO_barplot.pdf", width=7, height=7)
## Dot plot
p2 <- dotplot(ego, showCategory=10)
ggsave("GO_Dotplot.png", width=8, height=7)
ggsave("GO_Dotplot.pdf", width=8, height=7)
## Gene-Concept Network
p3 <- cnetplot(ego)
## Enrichment Map
#install.packages('ggnewscale')
library(enrichplot)
ego <- pairwise_termsim(ego)
p4 <- emapplot(ego)
分析结果展示
R版本信息
R version 4.1.2 (2021-11-01)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 22621)
2.2 topGO的安装和使用
与clusterprofiler一样，topGO也是一个发布在Bioconductor上的R包工具，与clusterprofiler不一样的是，topGO专注于GO富集分析，并且可以使用DAG图来可视化GO之间的结构关系
安装
if (!requireNamespace("topGO", quietly = TRUE))
BiocManager::install("topGO")
if (!requireNamespace("Rgraphviz", quietly = TRUE))
BiocManager::install("Rgraphviz")
数据准备
Gene2GOs.map - 目标物种的背景注释文件，总共两列，列与列中间用制表符（即Tab键）分隔，其中，每一个基因ID可能会对应多个GO ID (GO:xxxxxxx) ，格式如下：
gene_ID1 GO_ID1, GO_ID2, GO_ID3, ...
gene_ID2 GO_ID1, GO_ID2, GO_ID3, ...
gene_ID3 GO_ID1, GO_ID2, GO_ID3, ...
...
gene_ID100 GO_ID1, GO_ID2, GO_ID3, ...
GeneList.csv - 用于分析的目标基因文件，可选择只保留一列用于富集分析的基因ID，如果需要其余辅助信息，如fold change、pvalue等，最好将基因ID放在第一列，格式如下：
gene_ID1
gene_ID2
gene_ID3
...
gene_ID100
简单使用方法
# library packages
library(topGO)
library(Rgraphviz)
# 清除变量
rm(list=ls())
# 重新设置你的工作目录，你可以把它设置为分析数据文件所在的目录
setwd(".")
# step1: 读取数据
# 1 读取背景文件信息
geneID2GO <- readMappings(file = 'Gene2GOs.map')
geneNames <- names(geneID2GO)
# 2 读取目标基因文件
genes <- read.csv("GeneList.csv", header = FALSE, colClasses = c("character"))
myInterestingGenes <- genes[, 1]
geneList <- factor(as.integer(geneNames %in% myInterestingGenes))
names(geneList) <- geneNames
# 3 构建topGOdata对象
# 此处的`ontology`参数支持：BP，MF, CC
GOdata <- new("topGOdata", ontology = "BP", allGenes = geneList, annot = annFUN.gene2GO, gene2GO = geneID2GO)
# all GO terms are available for analysis
ug <- usedGO(GOdata)
# step2: 执行富集检验
# runTest functions return an object of type topGOresult
resultFis <- runTest(GOdata, algorithm = "classic", statistic = "fisher")
#pvalFis <- score(resultFis)
# step3: 结果分析和可视化
# 提取结果
allRes <- GenTable(GOdata, classic = resultFis, orderBy = "classic", ranksOf = "classic", topNodes = length(ug))
write.csv(allRes, paste("output", "topGO.csv", sep = "_"), row.names = FALSE)
# 可视化 GO structure
showSigOfNodes(GOdata, score(resultFis), firstSigNodes = 5, useInfo = 'all')
# 保存为PDF或PS格式文件
printGraph(GOdata, resultFis, firstSigNodes = 5, fn.prefix = "output", useInfo = "all", pdfSW = TRUE)
分析结果展示
R版本信息
R version 4.1.2 (2021-11-01)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 22621)
2.3 goatools的安装和使用
GOATOOLS是一个发布在Github上用于GO富集分析的Python工具，统计分析功能十分丰富，同时也可以绘制DAG图来可视化GO之间的结构关系
安装
使用pip安装
pip install goatools
使用conda安装
conda install -c bioconda goatools
数据准备
GO的结构关系和描述信息数据下载
wget http://current.geneontology.org/ontology/go-basic.obo
wget http://current.geneontology.org/ontology/subsets/goslim_generic.obo
association.txt - 目标物种的背景注释文件，总共两列，列与列中间用制表符（即Tab键）分隔，其中，每一个基因ID可能会对应多个GO ID (GO:xxxxxxx) ，格式如下：
gene_ID1 GO_ID1; GO_ID2; GO_ID3; ...
gene_ID2 GO_ID1; GO_ID2; GO_ID3; ...
gene_ID3 GO_ID1; GO_ID2; GO_ID3; ...
population.txt - 目标物种的背景基因文件，只有一列，格式如下：
gene_ID1
gene_ID2
gene_ID3
...
gene_IDN
study.txt - 用于分析的目标基因文件，比如100个差异基因，只保留一列，格式如下：
gene_ID1
gene_ID2
gene_ID3
...
gene_ID100
简单使用方法
# 执行GO富集分析
find_enrichment.py study.txt population.txt association.txt \
--pval=0.05 --method=fdr_bh --pval_field=fdr_bh \
--outfile=output_id2gos.xlsx,output_id2gos.tsv \
--obo=go-basic.obo \
--goslim=goslim_generic.obo \
--ns=BP,MF,CC \
--no_propagate_counts
# 绘制DAG图（有向无环图）
# 这里只选择了BP相关的GO id进行了绘制
sed '1d' output_id2gos.tsv | awk '$2=="BP"{print $1}' > BP_goids.txt
go_plot.py --go_file=BP_goids.txt --outfile=BP_goplot.pdf \
--id2gos=association.txt --obo=go-basic.obo --title=BP_DAG
分析结果展示
作者 ：LXD
转载 ：如何简单快速地完成一次GO富集分析
来源 ：微信公众号
著作权归作者所有，任何形式的转载都请联系作者。
