★ FEATURED ARTICLE 做植物转录组或基因组项目时第一步拿到序列后的固定动作就是功能注释。而在所有注释体系里GO注释几乎是必做的一环不管后续是要做差异基因富集、构建调控网络还是单纯在文章里放一张BarplotGO注释结果都是整个下游分析的“地基”。这篇就围绕植物基因组的GO注释把从原理、工具选型到命令实操、结果清洗再到富集可视化的完整流程拆开讲清楚全程以实际跑过的项目为例适合刚入门的生信同学也给正在优化注释流程的老手提供一些参数和坑位参考。1. GO注释是什么植物场景里为什么绕不开它1.1 GO的三层世界从分子功能到生物学过程GO的全称是Gene Ontology基因本体论它本质上是一套标准化的功能词典用来描述基因和基因产物的功能属性。它不是某一家数据库搞的私有分类而是由Gene Ontology Consortium维护的、跨物种统一的受控词汇表。这套词典把基因功能拆成三个完全不同的维度细胞组分Cellular ComponentCC基因产物在细胞里待的位置比如叶绿体类囊体膜、线粒体内膜、细胞壁。分子功能Molecular FunctionMF基因产物在分子层面的“手艺活”比如激酶活性、DNA结合、氧化还原酶活性。生物学过程Biological ProcessBP多个分子功能协作后完成的宏观过程比如光合作用、干旱响应、开花时间调控。为什么要分三层因为一个基因往往不只干一件事。比如一个典型的植物蛋白激酶它在细胞质里定位CC具有ATP结合和磷酸转移酶活性MF参与MAPK级联信号传导和防御反应BP。如果只给一个标签信息损失太大三个维度分开定义才能比较完整地描述基因的多面角色。1.2 植物GO注释的特殊之处别把动物流程直接搬过来很多做动物研究的同学转到植物上第一个不适应的就是GO注释的覆盖率和证据来源差异。植物基因组里存在大量动物基因组没有的基因家族比如NBS-LRR抗病基因、C4途径相关酶、苯丙烷代谢通路相关基因、光周期调控相关基因等。这些基因很多是植物谱系特有的通用数据库里的同源证据往往不足。再加上植物基因组本身的特点很多作物基因组大、重复序列比例高、基因复制事件频繁导致多基因家族成员之间序列相似性极高。比如小麦的六倍体基因组里同一个NLR基因可能有几十个拷贝这些拷贝的GO注释如果不能区分功能分化下游富集分析很容易出现冗余和偏差。所以植物GO注释的策略不能照搬人类或模式动物那套“调个blast2go就行”的思路需要在证据来源、参数阈值、物种特异性数据库几个层面做调整。这也是后面要重点展开的部分。1.3 拿到GO注释之后能干什么GO注释不是终点它是几乎所有下游功能分析的前提。最常见的用途有三个差异表达基因的功能概览拿到几百上千个差异基因先看它们整体偏向哪些BP、MF、CC形成初步生物学故事。GO富集分析统计哪些功能类别在目标基因集里显著富集这是文章里最常用的功能解读手段。构建功能网络和基因集打分比如GSEA、GSVA这类算法都需要预先对每个基因打上功能标签再结合表达量做集合水平的统计分析。一句话总结没有GO注释你手里永远只有一堆序列号什么生物学结论都讲不出来。2. 工具选型思路七个常用方案怎么挑GO注释的工具有很多但针对植物数据真正值得反复使用的就那么几个。下面把主流的方案拉一张表对比下再逐一说明各自的适用场景和坑。2.1 工具横向对比工具/数据库原理速度植物适用性上手难度典型产出InterProScan蛋白结构域/特征搜索较慢本地很好domain证据独立于同源中TSV带GO、InterPro条目eggNOG-mapper预计算直系同源簇映射快好支持数百个物种植物谱系较全低emapper.annotations带GO条目Blast2GO / OmicsBoxBLAST 数据库映射 注释规则极慢依赖BLAST好老牌方案可定制规则中高注释txt、GAF、富集图表AgriGO基于已知GO注释做富集也可做简单注释在线取决于网络专为植物设计低富集结果tableTBtools的GO注释模块封装eggNOG/InterPro等取决于底层引擎很好图形界面友好极低注释表格、富集结果PANNZER基于SVM和整合证据在线/本地中等一般微生物更强中GO TSVTrinotate基于转录本序列注释流程中可用于植物转录组非专门中SQLite数据库带GO2.2 实战中我的选型逻辑如果手头拿到的是拟南芥、水稻、玉米这类参考基因组注释蛋白我的首选是eggNOG-mapper直接跑速度快一次能出GO、KEGG、COG三套餐后面富集分析和通路注释都有了。如果是对一个全新测序基因组做从头注释得到的蛋白集我会用InterProScan跑一轮domain搜索作为主要证据因为InterProScan不依赖与已知物种的同源信息而是基于结构域模型对孤儿基因和新基因家族更友好。Blast2GO这几年我用得越来越少主要问题是它依赖BLAST步骤速度太慢而且受NCBI NR数据库的版本影响大。但如果项目要求手动精细调控注释规则比如设置严格的相似度阈值、EVALUE阈值Blast2GO仍然是最可控的方案。OmicsBox是它的商业升级版界面和注释后分析做得好只是授权费用不低很多课题组未必愿意掏。TBtools在植物领域其实已经成了“人手一个”的工具箱它的GO注释模块底层封装了InterProScan或eggNOG-mapper对完全不想碰命令行的人很友好。我个人建议是命令行方案还是要会因为服务器上没有图形界面脚本化的流程可以批量重跑参数试验这是GUI替代不了的。2.3 版本和数据库更新的隐患用工具之前一定先看版本。InterProScan每隔几个月就会更新数据库旧版本注释出的结果和新版本可能有较大差异特别是在结构域家族更新换代时会直接影响基因的GO term变化。eggNOG-mapper同样如此它有多个数据库版本如5.0、2024版等旧版本的物种覆盖和新版本差别不小新的植物基因组往往在更新版本里才有收录。我的习惯是每个项目的分析环境用conda固定版本同时记录下数据库的下载日期和版本号。写文章或报告时也建议注明版本信息否则别人无法复现你的结果审稿人问到也说不清楚。3. 实操流程从蛋白序列到GO注释结果3.1 准备输入文件不管用什么工具输入的常规格式都是FASTA蛋白序列不是核酸序列。常见错误就是直接把CDS序列丢进去这会导致注释结果里全是“翻译后修饰”或者“核糖体相关”这种模糊条目因为注释工具读的是氨基酸序列的保守特征。如果你手里只有CDS先用序列翻译工具转成蛋白序列。这一步要注意阅读框和终止密码子推荐用seqkit的翻译命令或者getorf等工具。实际操作# 用seqkit翻译CDS为蛋白序列剔除不完全的序列 seqkit translate cds.fa --frame 1 --trim --clean protein.fa如果是从基因注释GFF提取的CDS务必要保证CDS能完整翻译关于相位问题建议用gffread来处理更为稳妥。3.2 路线一eggNOG-mapper跑快速注释如果你选择eggNOG-mapper步骤非常简单conda装好后几乎只有两步。先下载数据库这一步很大通常需要几十GB空间建议留足磁盘空间再运行注释。# 下载eggNOG数据库新版可能需要先创建目录 mkdir -p eggnog_db cd eggnog_db download_eggnog_data.py -y --data_dir ./ # 运行注释注意输入格式是fasta emapper.py \ -i protein.fa \ --output prefix_plant \ --output_dir ./emapper_result \ --data_dir ./eggnog_db \ -m diamond \ --override \ --cpu 16关键参数说明-m diamond表示用DIAMOND做序列搜索这是eggNOG-mapper 2.0之后的主流方式比原版BLAST快两个数量级以上。--cpu按服务器核数调整理论上这个工具做序列相似性搜索时可以压缩到很小时间实际跑一个3万基因的植物蛋白集16线程大概几十分钟到两小时。运行结束后会生成多个文件核心是.emapper.annotations。这个文件里包含了每个基因对应的GO term、KEGG通路、COG分类号等。GO相关的主要是GOs这一列多个GO term用逗号分隔。值得注意的是有的eggNOG-mapper版本输出里GO term前缀会带GO:有的不带清洗时要注意统一。3.3 路线二InterProScan的结构域证据注释另一个更“硬核”的方案是InterProScan它搜索的不是同源序列而是蛋白结构域和功能位点数据库包括Pfam、PANTHER、SUPERFAMILY、Gene3D等等。对植物来说InterProScan和电子注释结果互补性强。# 安装可以直接通过conda/apt注意Java版本兼容 interproscan.sh \ -i protein.fa \ -f TSV \ -goterms \ -iprlookup \ -pa \ -o interpro_result.tsv \ -cpu 16参数解释-goterms让InterProScan解析结构域对应的GO条目-iprlookup在结果中附带InterPro条目ID-pa表示分析所有匹配对每个序列保留所有有意义的结构域而不是只保留best hit。这一步很重要因为植物蛋白往往由多个domain组成只留best hit会丢失大量信息。InterProScan的TSV输出每一行代表一个“序列-结构域-GO”的匹配关系一个基因可能出现多行。后面做基因到功能注释的去重时需要先按基因分组再把所有GO term合并起来去重。3.4 两条路线怎么合并实际项目中最稳的做法是同时跑eggNOG-mapper和InterProScan然后把两者的GO注释合并。因为eggNOG-mapper的同源映射速度快但依赖数据库覆盖InterProScan的结构域注释不依赖同源但对某些长非编码蛋白或新型蛋白可能漏判。两者合并后覆盖率通常能显著提高。合并时要注意去重和冲突处理同一个基因在两个来源都注释到同一个GO term取并集即可。不同来源给出不同GO term如果都在同一层级且不矛盾保留全部。如果一个来源注释到“分子功能未知”比如GO:0003674是MF的根节点另一个有具体条目时丢弃根节点条目。这一步可以写一个简短的Python脚本处理import csv from collections import defaultdict gene2go defaultdict(set) for source_file in [emapper_no_unknown.txt, interpro_go.txt]: with open(source_file) as f: for line in f: parts line.strip().split(\t) if len(parts) 2: continue gene, go_terms parts[0], parts[1] for go in go_terms.split(,): go go.strip() if go and not go.endswith(0003674): # 丢根节点 gene2go[gene].add(go) with open(merged_gene2go.tsv, w) as out: for gene, gos in gene2go.items(): out.write(f{gene}\t{,.join(sorted(gos))}\n)当然这是极简版实际项目还需要过滤有可疑来源的向上归类条目具体看工具注释规则后面会展开讲。4. 注释结果的清洗、评估与格式转换4.1 从注释表到gene2go三层拆分拿到合并结果后第一件事就是把所有GO term拆到BP、MF、CC三个维度下。因为GO的term是从根节点延伸下来的不区分维度去做富集分析会导致生物学含义混乱。最简单的方式是直接用go.obo文件解析每个GO term所属的namespace然后做映射。如果你不想写解析器可以直接用goatools这个Python库它内置了obo文件读取和GO term的层级关系查询。实际用法from goatools.obo_parser import GODag godag GODag(go.obo) with open(gene2go.tsv) as f, open(gene2go_split.tsv, w) as out: for line in f: gene, gos line.strip().split(\t) for go in gos.split(,): if go in godag: namespace godag[go].namespace out.write(f{gene}\t{go}\t{namespace}\n)输出文件里每一行一个“基因-GO-term-所属类”对应关系后续无论是用R语言的clusterProfiler做富集还是统计功能类别的分布都建议基于这个标准三元组表。4.2 覆盖率低的原因排查植物项目做完注释后最常见的焦虑就是“怎么这么多基因没注释到GO”。我遇到过拟南芥三万个蛋白里走了标准流程后还有三成没有GO的情况。排查思路从上到下序列质量问题转录组从头组装后翻译出的蛋白很多是片段结构域被截断同源搜索自然找不到。优先用TransDecoder或者保留ORF完整的序列。数据库版本太旧有些公共数据库里收录的植物物种有限如果你的物种比较小众比如苔藓、藻类、蕨类用旧版eggNOG可能覆盖很差。阈值设置不合理eggNOG-mapper支持调整--evalue和--score阈值默认值通常比较严格。可以放宽一些试试但要注意放宽后会带来假阳性。植物特有基因确实多很多非模式植物有大量的物种特有基因没有同源证据也没有domain信息这部分是客观盲区不用强求。4.3 统计报告怎么向老板或审稿人展示注释流程完成后一定要出一份注释率统计表。常见做法总基因数、有GO注释的基因数、注释率。分别统计BP、MF、CC三类的注释覆盖率。统计每个基因的GO term数量分布比如0个、1个、2-5个、6-10个、10个的基因数分别有多少。这类统计既能说明注释质量也能发现问题。比如如果一个基因平均GO term数是2以下说明可能存在上游注释偏严的问题如果注释率超过90%但大部分基因只有“binding”这类泛泛的MF条目说明证据质量不高下游分析一样会受影响。5. 下游进阶GO富集分析和结果解读5.1 富集分析的三种常见算法拿到gene2go表之后下一步通常是做GO富集分析。常用的方法有三类超几何分布/费希尔精确检验最简单适合比较一个目标基因集如差异基因和背景基因集如所有注释基因的差异显著性。代表工具是clusterProfiler的enrichGO。基于拓扑结构的富集如topGO的elim和weight算法考虑GO层级结构避免祖先节点冗余适合做系统化的层级分析。GSEA类基因集富集不设硬性差异阈值使用全部基因的表达值排序适合转录组中无显著差异但趋势一致的功能类别。对植物项目我习惯先跑clusterProfiler的enrichGO再用topGO做一次交叉验证。原因是两者对GO层级结构的处理方式不同如果两者都显著的结果通常非常稳如果只有某个算法报显著往往说明结果对算法敏感要慎重解读。5.2 用clusterProfiler跑植物GO富集R语言走起模拟一个水稻差异基因分析场景library(clusterProfiler) library(org.Osativa.eg.db) # 水稻注释包其他植物可用org.At.tair.db等 # 读取差异基因 deg - readLines(deg_list.txt) # 读取背景基因通常是所有检测到的基因 universe - readLines(background_genes.txt) # 读取gene2go映射表也可以用AnnotationDbi从注释包获取 gene2go - read.delim(gene2go.tsv, header FALSE) colnames(gene2go) - c(gene, go, namespace) # GO富集 ego_bp - enricher( gene deg, universe universe, TERM2GENE gene2go[, c(go, gene)], TERM2NAME go2term_table, # go term名称映射 pAdjustMethod BH, pvalueCutoff 0.05, qvalueCutoff 0.05 ) head(as.data.frame(ego_bp))这里最重要的一个前提是TERM2GENE必须是你自己的物种注释结果不要偷懒用模式物种的注释包。因为模式物种注释包里的基因ID和你项目里使用的ID不一定能对应上强行对应会丢大量基因富集结果完全失真。5.3 可视化从气泡图到网络图富集结果的可视化套路基本固定气泡图是文章标配横轴是GeneRatio或富集因子纵轴是GO term点大小代表基因数点颜色代表P值或Q值。网络图适合展示相互关联的GO term用emapplot可以直观看到哪些功能模块同时受到扰动。Reactome或Pathview这类通路图在植物里不常用因为植物通路注释远不如人类完善专心做GO富集的Bar图和气泡图已经能讲清楚大部分故事。一个容易忽视的细节作图时GO term名称应尽量展示“简洁可读”的短名而不是一串形如GO:0042742的ID否则读者根本不知道在讲什么。短名可以从go.obo提取def或name列也可以在clusterProfiler的结果表里直接取Description列。6. 植物GO注释的进阶策略与避坑心得6.1 多物种项目如何统一注释基准如果你做的是比较转录组牵涉到多个物种比如小麦族内部小麦、大麦、黑麦各物种的基因ID体系完全不一样但GO注释体系是统一的。此时建议不要分物种单独跑eggNOG-mapper后直接比较GO term因为各物种数据库覆盖度不一致会导致系统偏差。推荐方案是把所有物种的蛋白序列合并到一个FASTA里统一跑一次eggNOG-mapper或InterProScan保证每个物种的注释证据标准完全一致。后续在对比时再按物种拆开。这个思路能极大减少“拟南芥注释好、黑麦注释差”这类系统误差。6.2 GAF格式与公共数据库上传如果你打算把注释结果上传到公共数据库比如EBI的GOA或者和别人交换数据就必须把gene2go表转换成GAFGene Association File格式。GAF格式有严格的规定包括15列以上各行依次是数据库来源、基因标识、证据代码、参考来源等。不要自己硬造推荐使用goatools提供的GafWriter或者参考GOA的帮助文档生成。平时做项目可以不用GAF但如果你开发的是新物种参考基因组上传一套标准GAF格式的GO注释会给后续所有用这个基因组的人带来极大方便也算是对社区的一种回馈。6.3 证据代码被忽视但很重要的信息GO注释里还有一个证据代码字段Evidence Code比如IEAInferred from Electronic Annotation、ISSInferred from Sequence Similarity、IMPInferred from Mutant Phenotype。在植物自动注释流程中绝大部分基因都是IEA这是计算预测不是实验验证。写文章时务必备注清楚证据来源审稿人如果较真看到你拿IEA的注释结果去支撑一个强功能结论往往会质疑。如果有实验数据比如突变体表型、ChIP-seq靶基因可以单独整理一个“实验证据注释”子集标注成IDA/IMP/IPI等作为机器注释的补充和校正。6.4 时间成本与计算资源预估最后聊一聊资源和时间eggNOG-mapper 3万基因蛋白集16核推荐2~4小时以内能完成搜索和注释。InterProScan跑同样的数据会更耗时4~8小时都很正常取决于是否启用了全部数据库以及是否用并行的Precompute等。合并清洗、富集分析全流程下来视服务器负载和使用R/Python熟练度通常在1~2天内能全部跑完。内存方面eggNOG-mapper峰值占用不算高通常16~32GB足够InterProScan如果同时跑很多线程建议内存给到64GB以上否则容易进程被杀。数据库下载时注意网络稳定性最好挂后台下载或者用支持断点续传的工具。6.5 引文与报告规范做完一个完整的注释流程我记得要在文章的方法学部分写清楚以下内容GO数据库版本如go.obo是从哪天下载的版本号多少注释工具及版本号eggNOG-mapper 2.1.12, InterProScan 5.66-95.0等序列比对参数evalue阈值、score阈值、数据库版本物种分类与背景基因集定义富集分析的算法、FDR阈值、背景基因集合这五条写清楚别人基本能完整复现你的流程。很多文章直接写“We performed GO annotation using eggNOG-mapper”没有任何版本和参数信息这种做法在公共数据越来越受重视的今天已经不太够用了。7. 现场排查五个高频异常和处理办法实际跑植物GO注释时遇到的坑多到可以单独写一篇。这里挑五个高频的直接上排查结论。7.1 运行eggNOG-mapper提示数据库不存在新版eggNOG-mapper对数据目录要求严格--data_dir指向的目录需要有完整目录结构而且某些版本下载数据库后还会生成一个eggnog.db和一个eggnog_proteins.dmnd两者缺一不可。遇到数据库找不到优先检查这两个文件是否存在以及文件名大小写是否匹配。更常见的是并行下载时数据库文件损坏。建议用官方提供的校验命令或者在下载完成后跑一次--test自测。7.2 InterProScan报Java版本错误InterProScan 5以上依赖Java 11或17如果用系统自带的旧版Java启动即报错。解决方案就是显式指定JAVA_HOMEexport JAVA_HOME/path/to/jdk-17 export PATH$JAVA_HOME/bin:$PATH另外InterProScan解压包自带的interproscan.properties里有些路径不能有中文或空格放在服务器上时尽量用纯英文路径。7.3 注释率突然比别人的文章低一半不少人问我为什么自己的注释率只有50%人家文章里能到80%。先看物种如果做的是没有参考注释的非模式植物但对方是拟南芥、水稻那完全没法比。再看输入对方可能用的是UniProt的注释好的蛋白集而你用的是从头预测的基因模型两者标注完整度天然有差距。还有一个细节看对方是只统计“有任一GO注释”的基因还是要求三个namespace都至少有一个条目不同口径差异很大写文章和对比时务必统一统计口径。7.4 GO term里出现大量“intracellular”和“cellular process”这类过于宽泛的条目虽然技术正确但信息量接近于零在富集分析里还容易被误当成显著项。处理方式通常是基于GO层次结构设置最低深度depth过滤比如只保留从根节点往下至少三层或四层的条目。用goatools或者GO层次结构数据能做到。这个处理和“IHOP”无关纯粹是注释质量控制的实操技巧。我在项目里会写一个过滤脚本把每个GO term的注释深度、子节点数量都算出来统一做一次去泛化。7.5 富集结果全是三羧酸循环和翻译没有植物抗逆相关条目排除样本问题后这种结果大概率是背景基因集污染了。比如你拿差异基因做富集但背景基因里混入了很多从NCBI下载的、与当前物种无关的基因ID导致背景功能分布不均匀富集信号被误导。解决方法是严格定义背景为“本次实验检测到表达的所有基因”而不是“物种全部基因”更不是“注释数据库全部基因”。8. 写在最后顺手可用的几个改进思路GO注释这件事表面看就是个跑工具的过程但真正决定注释质量的是对物种特征、工具原理和数据质量的把握。这里再说几个我自己实践后觉得有用的思路。第一尽量把GO注释和KEGG注释放在同一个流程里跑。eggNOG-mapper同时会产出KEGG通路注释一次搜索同时拿到两套结果后续富集分析时还能做GO-KEGG的交叉引用节省大量重复计算时间。第二针对植物特异性数据库做一些定制。比如做茄科作物可以考虑额外补充SGN数据库的注释信息做禾本科可以补充Gramene或Ensembl Plants的已有注释不必完全依赖默认数据库。这些物种特异性补充往往能明显改善稀有基因家族的注释覆盖率。第三不要忽略“无注释”基因的价值。有些基因没有GO注释不代表它们没有功能。在植物基因组里新基因、快速演化的物种特异基因在GO注释里经常是缺失的反而这些基因可能是物种适应性进化的关键。可以单独对这些基因做保守性分析、结构域扫描和非同义突变率分析说不定能挖出新的研究亮点。我在植物生信项目里跑了这么多轮GO注释最深的一个体会是注释不是“一键生成结果”的事每一个GO term背后都有证据来源、有层级关系、有阈值影响拿到结果后花时间做质控比盲目追求更高的注释率更值得。希望你在这份流程基础上能跑出一份经得起审稿人和合作者检验的注释结果。
阅读完成 · 觉得有帮助？
