基因本体GO注释详解
一句话概述：Gene Ontology（GO）是描述基因产物功能的标准化词汇体系，分为生物过程(BP)、分子功能(MF)、细胞组分(CC)三大类，是功能注释和富集分析的基础。
核心知识点速查表
| 概念 | 说明 |
|---|
| GO | Gene Ontology，基因本体论（白话：给基因功能贴标签的标准体系） |
| BP(Biological Process) | 生物过程（白话：基因参与了什么活动，如"细胞分裂"） |
| MF(Molecular Function) | 分子功能（白话：基因产物能干什么，如"激酶活性"） |
| CC(Cellular Component) | 细胞组分（白话：基因产物在细胞哪个位置，如"线粒体"） |
| GO term | 具体的功能描述条目，如GO:0006915(apoptotic process) |
| DAG | 有向无环图，GO的层级结构（从泛到精越来越具体） |
| ORA | Over-Representation Analysis，超几何检验富集分析 |
| GSEA | Gene Set Enrichment Analysis，基因集富集分析 |
| FDR | False Discovery Rate，多重检验校正后的p值 |
一、GO注释方法
1.1 使用InterProScan注释
# === InterProScan：最全面的GO注释工具 ===
# 整合了Pfam、SMART、CDD、PANTHER等多个数据库
# 安装
conda install -c bioconda interproscan # conda安装
# 运行
interproscan.sh \
-i proteins.fasta \ # 输入蛋白质序列
-f tsv,gff3 \ # 输出格式（表格+GFF3）
-goterms \ # 输出GO注释
-pa \ # 输出通路注释
-dp \ # 禁用前查看（直接运行）
-cpu 16 \ # CPU线程数
-o interpro_results # 输出前缀
1.2 使用eggNOG-mapper注释
# === eggNOG-mapper：快速功能注释 ===
# 基于预计算的直系同源组，速度快
# 安装
conda install -c bioconda eggnog-mapper
# 下载数据库
download_eggnog_data.py --data_dir eggnog_db # 下载数据库（约50GB）
# 运行
emapper.py \
-i proteins.fasta \ # 输入蛋白质序列
--data_dir eggnog_db/ \ # 数据库目录
--output eggnog_out \ # 输出前缀
--cpu 16 \ # CPU线程
-m diamond \ # 使用DIAMOND搜索（更快）
--go_evidence all \ # 输出所有证据水平的GO
--tax_scope auto # 自动选择分类范围
1.3 使用Blast2GO（商业工具替代方案）
# === 用BLAST + Mapping方式获取GO注释 ===
# 第1步：BLAST搜索UniProt/SwissProt
blastp \
-query proteins.fasta \ # 输入蛋白质
-db uniprot_sprot \ # SwissProt数据库
-outfmt "6 qseqid sseqid pident evalue" \ # 表格输出
-evalue 1e-5 \ # E-value阈值
-max_target_seqs 5 \ # 每个query保留5个hit
-num_threads 16 \ # 线程数
-out blast_results.txt # 输出文件
# 第2步：从UniProt ID映射GO（Python）
python3 << 'PYEOF'
import requests
import pandas as pd
# 读取BLAST结果
blast = pd.read_csv("blast_results.txt", sep="\t",
names=["query", "subject", "identity", "evalue"])
# 提取UniProt ID
uniprot_ids = blast["subject"].str.extract(r'sp\|(\w+)\|')[0].unique()
# 批量查询UniProt获取GO
for uid in uniprot_ids[:10]: # 示例：查前10个
url = f"https://rest.uniprot.org/uniprotkb/{uid}?fields=go_id,go"
resp = requests.get(url)
if resp.ok:
print(f"{uid}: {resp.json().get('go', 'No GO')}")
PYEOF
二、GO富集分析
2.1 R语言clusterProfiler
# === 使用clusterProfiler做GO富集分析 ===
library(clusterProfiler) # 加载包
library(org.Hs.eg.db) # 人类注释数据库（换成对应物种）
# 准备基因列表
gene_list <- read.table("diff_genes.txt", header=FALSE)$V1 # 差异基因ID
# ORA超几何检验富集
ego <- enrichGO(
gene = gene_list, # 差异基因列表
OrgDb = org.Hs.eg.db, # 物种注释数据库
keyType = "ENTREZID", # 基因ID类型
ont = "ALL", # BP/MF/CC/ALL
pAdjustMethod = "BH", # 多重检验校正方法
pvalueCutoff = 0.05, # p值阈值
qvalueCutoff = 0.05 # FDR阈值
)
# 查看结果
head(as.data.frame(ego), 20) # 前20条富集结果
# 可视化
dotplot(ego, showCategory=20) # 点图（最常用）
barplot(ego, showCategory=15) # 柱状图
cnetplot(ego, showCategory=5) # 网络图
goplot(ego) # GO DAG图
2.2 GSEA分析
# === GSEA：不需要设阈值的富集分析 ===
library(clusterProfiler)
library(org.Hs.eg.db)
# 准备排序的基因列表（按logFC或统计量排序）
gene_rank <- read.table("gene_rank.txt", header=TRUE) # 基因ID + 排序值
geneList <- gene_rank$logFC # 排序值
names(geneList) <- gene_rank$ENTREZID # 基因ID
geneList <- sort(geneList, decreasing=TRUE) # 从大到小排序
# 运行GSEA
gsea_go <- gseGO(
geneList = geneList, # 排序的基因列表
OrgDb = org.Hs.eg.db, # 注释数据库
keyType = "ENTREZID", # ID类型
ont = "BP", # 只看生物过程
pAdjustMethod = "BH", # FDR校正
pvalueCutoff = 0.05 # 阈值
)
# 可视化
gseaplot2(gsea_go, geneSetID=1:3) # GSEA曲线图（前3个通路）
ridgeplot(gsea_go, showCategory=15) # 山脊图
2.3 非模式生物GO富集
# === 非模式生物：自定义GO注释 ===
# 用eggNOG-mapper或InterProScan的结果
# 第1步：准备gene2go映射文件
# 格式: GeneID GO_term
# gene001 GO:0006915
# gene001 GO:0005737
# gene002 GO:0006468
# 第2步：读取映射
gene2go <- read.table("gene2go.txt", header=FALSE,
col.names=c("gene", "GO"))
# 构建term2gene格式
term2gene <- gene2go[, c("GO", "gene")] # 第一列GO, 第二列gene
# 第3步：获取GO term描述
library(GO.db)
term2name <- select(GO.db, keys=unique(term2gene$GO),
columns="TERM", keytype="GOID")
colnames(term2name) <- c("GO", "name")
# 第4步：富集分析
ego_custom <- enricher(
gene = diff_genes, # 差异基因列表
TERM2GENE = term2gene, # GO-基因映射
TERM2NAME = term2name, # GO ID-名称映射
pvalueCutoff = 0.05 # p值阈值
)
dotplot(ego_custom, showCategory=20) # 可视化
三、GO分析结果解读
关键指标
| 列名 | 说明 |
|---|
| Description | GO term描述 |
| GeneRatio | 差异基因中注释到该GO的比例 |
| BgRatio | 背景基因中注释到该GO的比例 |
| pvalue | 原始p值 |
| p.adjust | FDR校正后的p值（看这个！） |
| Count | 注释到该GO的差异基因数 |
| geneID | 具体的基因列表 |
解读原则
- 看p.adjust而非pvalue：多重检验后<0.05才可靠
- 关注GeneRatio：富集倍数高才有生物学意义
- BP最重要：通常先看生物过程，再看分子功能和细胞组分
- 结合背景知识：GO富集结果需要结合具体生物学问题解读
四、面试高频考点
Q1: GO的三个分类是什么？各举一例
- BP(Biological Process)：如GO:0006915 细胞凋亡过程
- MF(Molecular Function)：如GO:0005524 ATP结合
- CC(Cellular Component)：如GO:0005739 线粒体
- 白话：BP=参与什么活动，MF=有什么技能，CC=在哪工作
Q2: ORA和GSEA的区别？
| ORA | GSEA |
|---|
| 输入 | 差异基因列表（需要阈值） | 所有基因的排序列表 |
| 方法 | 超几何检验 | 排列检验 |
| 优点 | 简单直观 | 不依赖阈值，信息利用更充分 |
| 缺点 | 依赖阈值选择 | 计算量大，结果解释复杂 |
Q3: GO富集分析常见陷阱？
- 背景基因选错：应该用所有检测到的基因，不是全基因组
- 忽略多重检验校正：用p.adjust而非pvalue
- 过度解读通用GO：如"protein binding"太泛，无实际意义
- 忽略层级结构：父子GO term可能高度冗余
常见报错与解决
| 报错 | 原因 | 解决方案 |
|---|
No gene can be mapped | 基因ID类型不匹配 | 用bitr()转换ID类型 |
No enrichment found | 差异基因太少或注释不全 | 放宽阈值或用自定义注释 |
OrgDb not available | 非模式生物无OrgDb | 用enricher()自定义分析 |
Too many terms | 结果太多 | 使用simplify()去冗余 |
速查表
# === GO分析速查（R语言） ===
library(clusterProfiler)
library(org.Hs.eg.db) # 换成对应物种
# ORA富集
ego <- enrichGO(gene=gene_list, OrgDb=org.Hs.eg.db,
keyType="ENTREZID", ont="ALL")
dotplot(ego, showCategory=20)
# GSEA
gsea <- gseGO(geneList=sorted_genes, OrgDb=org.Hs.eg.db, ont="BP")
gseaplot2(gsea, geneSetID=1:3)
# 非模式生物
ego <- enricher(gene=genes, TERM2GENE=term2gene, TERM2NAME=term2name)
# ID转换
bitr(genes, fromType="SYMBOL", toType="ENTREZID", OrgDb=org.Hs.eg.db)
# 去冗余
simplify(ego, cutoff=0.7, by="p.adjust")