上次帮一个做环境微生物的朋友注释一批宏基因组bin他图省事直接拿真核注释流程去跑结果那叫一个惨——基因结构预测全乱套tRNA还缺了一大堆。后来换成prokka同样的数据十分钟跑完产物干净得可以直接往下游分析。prokka这个工具做微生物组或者细菌基因组研究的人应该都听过但真正用顺手的人不多。它是英国Wellcome Sanger Institute的Torsten Seemann开发的专门用来做原核生物和病毒的基因组自动注释。所谓注释说白了就是拿到一段FASTA格式的基因组序列之后告诉你有多少个基因、每个基因是干嘛的、编码什么蛋白、有哪些RNA元件。听起来简单但背后涉及基因预测、同源比对、功能注释好几层工作prokka的价值就是把这些步骤全部串成一条流水线一条命令跑完。这篇文章我不打算念手册就按我自己实际使用的经验把prokka从选型、安装、参数调到结果验证讲一遍中间会穿插我踩过的坑和现在固定的操作习惯。不管你是刚入门想给手头的细菌基因组做第一个注释还是已经跑过几次但总觉得结果哪里不太对这篇都值得看一下。1. 选型之前先搞清楚prokka到底在做什么什么场景别用它1.1 注释的本质把“天书”翻译成有生物学意义的语言一条细菌基因组序列在电脑里就是ATCG四个字母排成的长串长度通常在2Mb到10Mb之间。直接看这段序列你什么信息都得不到。注释的过程就是用算法在序列上找特征——哪里是基因的起始和终止位置基因之间有没有间隔哪些区域编码tRNA、rRNA然后把这些特征跟已知功能的数据库比对得出“这个基因可能编码某一种酶”之类的结论。prokka做的就是这个翻译工作它的核心流程分四步第一步用Prodigal做基因预测把可能的编码序列CDS找出来第二步分别用Aragorn和Barrnap预测tRNA和rRNA第三步把预测出来的蛋白序列跟公共数据库比对寻找同源信息第四步综合所有信息生成标准格式的注释文件。这套流程听起来不复杂但每一步都有不少细节坑。prokka的优势在于把整个流程做了封装参数选项合理输出格式规范还有个很实用的设计——它的输出文件会按GFF3、GenBank、FASTA等多种格式同时生成直接把下游分析的接口留好了。1.2 与主流工具对比prokka的生态位在哪里我见过不少人在选注释工具时纠结这里直接跟大家说结论现在原核注释领域主流就是prokka、RAST、PGAP和bakta这几个各有各的定位。RAST走得早以前是网页版为主现在也提供API和本地版但它的问题是比较老派注释结果里经常会混进一些来自模式菌株的错误基因模型而且数据库更新慢。PGAP是NCBI官方的pipeline注释质量理论上最高但它要求输入数据干净且无法对病毒和质粒做灵活调整运行速度也慢全流程加人工审核常常需要几天。bakta是后起之秀数据库更新快运行也快在某些情况下比prokka的预测准确率还高不过它对输入序列的完整度有要求contig太碎的时候表现一般。prokka的最大竞争力在一个“快”字上——单条细菌基因组16线程下一般20分钟以内能跑完而且结果稳定、格式标准、易上手。还有个细节是很多下游工具比如泛基因组分析工具Roary原生的输入格式就是prokka的GFF文件生态衔接特别顺畅。需要提醒的是prokka只适用于原核生物和病毒。它是基于原核基因结构的统计特征做预测的用的是Prodigal和Barrnap这些专门针对细菌古菌设计的程序。如果你手里是植物、动物或者真菌的真核基因组比如最近很多人讨论的芍药基因组这类就别用prokka了建议去用MAKER、BRAKER或者Augustus EvidenceModeler这套流程。真核基因有内含子、可变剪接原核注释工具完全处理不了。1.3 什么情况下不建议用prokka除了真核基因组应该绕行之外还有两类情况我也不建议硬上prokka。第一类是基因组质量太差、contig数量特别多、N50特别低的数据。prokka设计理念是“注释整个基因组”如果一条基因组被拆成几百甚至上千个contig基因预测时会因为序列断裂而产生大量截短的假基因注释结果的生物学意义就大打折扣了。这种情况我会建议先做组装质量的评估和优化把N50拉到比较合理的水平再注释不要抱着“先注释看看”的心态浪费时间和磁盘。第二类是注释高度特化的菌株或者人工改造过基因组的样本。prokka本质上是靠数据库比对来给基因挂功能标签的如果样本里有大量novel gene或者经过大规模基因编辑的序列比对结果会以“hypothetical protein”居多注释率偏低。这类数据更适合用PGAP配合人工审核甚至需要结合RNA-seq数据进行手工校正。2. 环境准备与安装我试过的两种方法直接给你结论2.1 方法一conda安装推荐prokka的依赖比较多包括Prodigal、Barrnap、Aragorn、Infernal、Minced、BLAST、HMMER等十来个组件。自己一个个装非常痛苦还会遇到版本兼容问题。所以我强烈建议用conda一条命令解决所有依赖。# 创建独立环境避免依赖冲突 conda create -n prokka -c conda-forge -c bioconda prokka # 激活环境 conda activate prokka # 初始化数据库索引 prokka --setupdb这里要说明一下prokka从1.14版本开始就建议安装后先跑一次prokka --setupdb它的作用是生成BLAST数据库索引和HMMER模型的hmmpress格式文件。如果跳过这步运行时经常会在barrnap或者BLAST比对环节报错。我现在的习惯是额外装一个mamba来管理生物信息环境因为conda在解决某些依赖碰撞时会慢得让人崩溃。用mamba替代conda执行同样的命令速度会快非常多mamba create -n prokka -c conda-forge -c bioconda prokkaconda-forge和bioconda这两个channel的顺序建议按我上面写的来不要调换不然可能拉到旧版本包。2.2 方法二源码安装或直接拷贝适合服务器离线环境实验室服务器经常是离线状态或者因为安全策略不让随便用conda。这种情况下我试过两种方式。一种是从GitHub拉源码自己配依赖git clone https://github.com/tseemann/prokka.git cd prokka/bin # 将prokka所在路径加入环境变量 export PATH$PATH:/path/to/prokka/bin但这种方式你必须自己保证所有依赖工具在PATH里找得到prokka的安装检查器prokka --check会列出缺失的组件但是逐个解决仍然心累。另一种方式更省事找一台能联网的机器配好conda环境把整个环境tar打包拷到目标机器解压后改一下环境路径就能用。这在生物信息学领域很常见我自己在迁移到计算节点的时候就这么干过。不过前提是操作系统架构要一致CentOS 7打包的环境不能直接拷到Ubuntu 22.04上glibc版本差异会导致各种诡异报错。2.3 数据库文件配置默认库就够用吗prokka自带的数据库在安装目录下的db文件夹里包含Kingdoms、Bacteria、Viruses等子目录。需要注意的是prokka的默认数据库其实只是“种子”数据库它并没有包含RefSeq或UniProt的全量数据。这意味着它对注释一些公共数据库中较少见的物种时命中率不会那么高。系统提供了prokka-genbank_to_fasta_db和prokka-uniprot_to_fasta_db两个辅助脚本允许你把自己的参考数据库转换成prokka可用的格式。实操中我一般会额外加一个由近缘物种蛋白质序列构建的小型参考库作为--proteins参数传给prokka这样注释精确度会明显提升。具体的用法在后面参数部分会展开。3. 核心参数与完整注释流程从最简命令到工业级用法3.1 一条命令搞定基本注释先看最基础也是最常用的命令格式。假设你手里有一个细菌的完整基因组序列文件strainA.fasta需要做注释命令如下prokka strainA.fasta --outdir annotation --prefix strainA --kingdom Bacteria --cpus 8跑完之后在annotation目录下会生成多个文件其中最关键的是strainA.gffGFF3格式的注释文件包含所有基因、RNA、CDS的位置和属性描述。下游分析比如Roary泛基因组分析就用这个。strainA.gbkGenBank格式可以直接用SnapGene、Artemis等工具图形化查看。strainA.ffn所有基因的核苷酸序列。strainA.faa所有基因对应的蛋白序列。strainA.txt注释统计摘要包含CDS数量、tRNA数量、rRNA数量等。这一条命令对很多场景已经够用了。但如果你只满足于这一步prokka至少七成的能力没被发挥出来。下面我把几个我认为值得关注的参数展开讲。3.2 关键参数逐一拆解每个参数背后都有我踩过的坑--kingdom这个参数用来指定注释物种的分类学范围可选项是Bacteria、Archaea、Mitochondria、Viruses。默认值是Bacteria。如果注释对象是古菌比如产甲烷菌这类不指定Archaea可能会导致tRNA扫描程序Barrnap选错模型把一些真细菌特有的RNA特征套到古菌序列上结果就是tRNA数量离谱。同理注释噬菌体的时候要指定Viruses因为病毒基因组结构跟细菌差异非常大基因密度和编码策略都不同。有一个很容易忽略的使用细节——线粒体基因组也可以指定--kingdom Mitochondria。很多人不知道这一点真核基因组里的线粒体单独提出来用prokka注释效率比套MAKER流程高很多。不过严格来说线粒体基因组的注释最好配合--gcode 2把密码子表改成线粒体遗传密码表否则翻译出来的蛋白会有一堆看起来很奇怪的长度。--genus和--species这个参数除了在产物的“locustag”里体现之外还间接控制了基因命名的方式。比如说--genus Salmonella --species enterica的情况下输出的基因名会是SALMONEL_00001这样的格式而不是默认的STRAINA_00001。从实际使用的角度这个参数的意义在于让后续结果里有可读性更好的标签。我在做批量物种注释时会把属名提出来作为前缀结果文件再合并的时候一眼就能分辨每个基因来自哪个样本不需要额外做映射表。--proteins这是整个prokka里我最推荐大家花时间研究的参数它可以接收一个FASTA格式的蛋白质序列文件作为参考数据库供最终注释时比对。它有三层作用一是提升注释的准确率二是让更多的基因能被挂上功能描述名称三是让产物中的EC编号、Gene Ontology条目更丰富。实际操作中我一般用两种方式来构建参考蛋白文件# 方式一从NCBI下载近缘物种的蛋白序列 # 假设已经下载了refseq_protein.faa prokka strainA.fasta --proteins refseq_protein.faa --outdir annotation # 方式二增量式自身的prokka结果迭代 # 第一次先用默认库跑出结果把.faa作为参考再跑一次 prokka strainA.fasta --proteins first_pass.faa --outdir annotation第二种方式看着有点绕但实际上效果很好。第一轮跑完后得到初步的蛋白序列集合第二轮以这些蛋白序列作为参考prokka会给更多CDS挂上具体的功能名称而不是“hypothetical protein”。我做过一个测试第二轮注释后hypothetical protein比例能从30%降到15%左右效果惊人。这个“用自身结果迭代优化注释”的思路很多教程里都没写过算是比较小众的经验。--rfam这个参数会让prokka额外用Infernal程序去搜索Rfam数据库识别非编码RNAncRNA。Rfam数据库收录了大量核酶、核开关riboswitch、小调节RNA等非编码RNA家族默认的barrnap只负责rRNA和tRNA识别不到这些。代价是速度变慢。因为Infernal用的是profile-HMM加共线性模型的比对方法比普通的BLAST要慢不少。我实测过在一条5Mb的细菌基因组上加上--rfam之后运行时间大约会增加40%~60%。如果做的是比较标准的细菌基因组建议还是加上核糖开关在代谢通路研究中经常有重要角色。如果是想对大批量的contig快速注释做初筛可以不加先跑完再筛目标菌株做精细注释。--cpusprokka的多线程不是传统意义上的并行它内部会把任务拆分成若干个子任务部分子任务用BLAST比对时支持多线程但整体加速不是线性的。我试过8线程和16线程对单条基因组的加速比大约1.8倍并没有翻倍的收益因为流程里还包含不少单线程步骤。一个更实际的经验是多任务并行时不要把每个任务都指定14、16个线程。比如一台64核的服务器为了跑多个样本我通常会给每个prokka任务分配4~6个线程同时开10个任务总体吞吐量反而比一次跑一个用48线程的任务高得多。这个调度经验在做批量注释时省了我非常多的时间。--gcode这个参数指定遗传密码表编号默认是11细菌、古菌的标准密码子表。但如果你注释支原体、某些纤毛虫等使用替代密码子表的物种就必须指定正确的密码子表。支原体用密码子表4也就是UGA编码色氨酸而不是终止信号不指定的话Prodigal会把大量正常的末端外显子识别错误导致一堆截短蛋白。我们做微生物的遇到支原体的情况不少这个参数是一个很容易被忽略但影响巨大的调节项。3.3 从原始序列到注释报告的完整实操下面用我最近处理的一个例子把全流程串起来。假设我从测序公司拿到一个组装好的枯草芽孢杆菌基因组文件叫subtilis.fasta。首先我会先看一眼序列的基础信息# 统计序列条数和总长度 grep -c ^ subtilis.fasta # 通常完整细菌基因组应该是1条序列如果是contig形式数量会更多 # 检查序列是不是有非法字符 grep -v ^ subtilis.fasta | grep -iE [^ATCGN]这一步是为了防止后面跑到一半才发现输入文件有问题。然后执行注释prokka subtilis.fasta \ --outdir subtilis_annotation \ --prefix subtilis \ --kingdom Bacteria \ --genus Bacillus \ --species subtilis \ --rfam \ --cpus 8 \ --proteins bacillus_ref_proteins.faa这里额外用了--proteins参考蛋白选的是Bacillus属已知功能的蛋白集合。跑的过程中终端会打印各阶段的进度日志包括“Predicting genes”、“Running RNA searches”、“Functional annotation”等步骤。正常情况下等待15~30分钟就会看到Finished字样。输出目录结构如下subtilis_annotation/ ├── subtilis.ffn # 基因的核苷酸序列 ├── subtilis.fna # 输入基因组的副本拷贝 ├── subtilis.faa # 蛋白序列 ├── subtilis.gbk # GenBank格式 ├── subtilis.gff # GFF3格式下游主流输入 ├── subtilis.log # 运行日志报错排查第一现场 ├── subtilis.png # 基因组圈图预览 ├── subtilis.sqn # 用于提交NCBI的Sequin文件 ├── subtilis.tbl # 用于提交NCBI的第2个文件 └── subtilis.txt # 注释统计摘要其中subtilis.png这个文件很多人没注意到它是prokka用CGView脚本自动生成的一个基因组圈图预览看着不那么精致但在快速确认基因分布密度是否均匀、有没有大片空白区域时很有用。3.4 输出文件逐一详解哪些是核心产物哪些要被下游工具消费拿到输出文件夹时新手常见的问题是搞不清该用哪个文件。我来按用途把文件分组第一组结果浏览型用subtilis.gbk。这个文件可以直接拖进SnapGene或者Geneious里查看基因结构、阅读CDS的方向和注释信息。做PCR引物设计的时候我也习惯从这里面提取目标基因的上下游序列。第二组下游分析型用subtilis.gff。Roary泛基因组分析、Prokka注释结果转成Tree比如用GToTree构建进化树都认这个格式。有一点要提醒prokka输出的GFF3里FASTA序列部分是直接拼在文件尾部的如果下游工具解析严格可能报格式错误。遇到这种报错先别怀疑prokka试着用agat_sp_extract_sequences.pl或者简单脚本把序列部分剥掉再传入工具。第三组序列提取型用subtilis.ffn和subtilis.faa。前者适合做同源基因搜索、设计探针后者适合做蛋白功能预测、构建物种系统发育树。第四组数据提交型用subtilis.sqn和subtilis.tbl。如果你要把基因组提交到NCBI的GenBank数据库这两个文件可以直接配合表格编辑器使用。但提醒一句NCBI近年的提交流程主要走Submission Portalprokka生成的sqn文件得经过Tabl2asn或BankIt处理才能顺利进入正式流程直接裸传不一定被接受。4. 注释质量验证跑完不等于跑对这四步帮你把关4.1 第一步看统计数字快速判断异常打开subtilis.txt文件关注几个核心指标CDS数量、tRNA数量、rRNA数量、假基因数量、注释到功能的比例。以枯草芽孢杆菌为例一个完整的基因组通常有4100~4500个CDStRNA在70~90个左右rRNA基因通常有10个对应5S、16S、23S各多个拷贝。如果你的样本注释出来tRNA只有20个或者CDS数量比预期少了一两千那大概率是基因组组装不完整或者参数设置有问题需要回去检查输入数据。之前有朋友注释一个肠杆菌科的菌株跑出来CDS只有1800个我让他去检查N50结果发现原始数据被错误地切掉了大量片段。这类问题prokka本身不会报错所以注释完一定要自己核对统计数字。4.2 第二步随机抽查基因的功能注释从注释结果的subtilis.faa文件里随机挑若干个蛋白跑一次BLAST搜索看看它们的功能描述是否跟已知的相似物种蛋白对得上。注意如果在NCBI上比对显示的同源蛋白所属物种跟你的样本亲缘关系差得很远比如你注释的是变形菌门的菌却比对出了厚壁菌门的蛋白那就要警惕数据库污染或者序列样本本身有交叉污染的问题。我自己的习惯是单次注释至少抽查20个CDS特别关注那些带“multi-domain protein”或“hypothetical protein”标签的序列。后者多了不要慌很多时候是数据库覆盖不全不代表预测错误。4.3 第三步与近缘参考菌株的注释做比较如果同属里已经有近缘物种的公共注释结果可以做一个宏观对比。比较CDS平均长度、基因密度、GC含量分布这些指标在两个近缘物种之间应该非常接近。我遇到过一次比较极端的情况一个样本注释出来的CDS平均长度比参考菌株短了300bp后来查出来是组装过程中混入了载体序列测序接头没切干净。这类问题靠检查原始数据才能定位但注释结果的异常通常能给你提供最早的线索。4.4 第四步检查关键功能基因是否存在如果实验背景比较明确你大概知道样本里应该有某些功能特征比如抗性基因、特定代谢通路基因那就直接在注释结果里搜一下这些关键词。# 在gff文件里找特定注释关键词 grep -i tetracycline subtilis.gff grep -i cellulase subtilis.gff这个方法看起来朴素但非常有效。我曾经在处理土壤宏基因组的分箱结果时就靠这一步快速确认了某个bin是否含有产甲烷通路的关键酶节省了大量时间。5. 常见报错与排查技巧把我踩过的坑直接说给你5.1 高频报错速查表报错信息可能原因解决方案ERROR: Could not run barrnapBarrnap依赖的HMM模型未初始化重新运行prokka --setupdbERROR: Problem with GFF3 file输入的GFF3格式不规范检查文件是否被Excel等工具损坏用head确认格式ERROR: Unable to locate Prodigal依赖程序不在PATH中检查conda环境激活状态用which prodigal验证WARNING: Skipping contig with Nscontig含有大量未解碱基N不必处理一般不影响整体注释ERROR: Cannot open FASTA file文件路径错误或权限问题检查文件权限ls -l确认无特殊字符路径ERROR: Sequence names not uniqueFASTA头部有重复序列名用脚本给所有contig改名确保唯一性补充一个我遇到的比较少见的坑——如果prokka --setupdb时报权限错误往往是因为conda环境目录存在写权限问题。最简单的说服方法是用root身份重新执行一次或者把conda环境安装在用户有写权限的目录下。5.2 内存炸掉怎么办大基因组的prokka优化prokka本身比很多组装工具要省内存但处理超大基因组比如某些Streptomyces基因组可达10Mb以上的时候BLAST比对阶段内存会明显上涨我见过峰值超过16GB的情况。如果内存吃紧有两个优化方向。第一给--proteins提供高度冗余的近缘参考库这样命中率上升后候选序列转BLAST的阶段需要处理的序列数量会减少。第二把--evalue调高到1e-06甚至1e-05过滤掉更多弱匹配比对时间会缩短。但要注意evalue太宽松会引入假阳性注释质量会下降1e-06是我比较推荐的平衡点。5.3 让prokka更快三个立竿见影的技巧第一个技巧在前面提过就是批量样本时用4~6线程开多个任务而不是单任务吃满所有核。第二个技巧是分阶段注释——把基因组切成几个片段分别跑prokka再合并注释结果这个操作在遇到超大基因组时效果好但要注意跨切割边界基因的处理一般我不太推荐新手去做。第三个技巧是选择性地关闭--rfam在初筛阶段省时间只在最终精细注释时打开。5.4 干扰因素排查版本更新和数据库污染prokka迭代比较快隔几个月版本更新后注释结果可能会产生细微差异。如果你要拿注释结果发表或者做多批次对比强烈建议在论文方法部分写明prokka的版本号并把conda环境的配置文件conda env export一起留档。我在实际工作中就踩过这个坑——换了一台服务器后环境版本不同同样的基因组注释出的CDS数量差了几十个后来花了不少时间排查历史版本环境才搞定。关于数据库污染的问题prokka不像某些商业软件那样内置严格的数据库清洗流程如果你自己往--proteins里加了来源不明的序列很容易引入错误的注释。建议所有自建参考库都做一次CD-HIT去冗余去掉同源性过高的序列保持参考库的干净度。6. 从注释往下游走prokka结果的典型应用和生态衔接6.1 泛基因组分析Roary的黄金搭档prokka跟Roary之间几乎是“无缝衔接”的关系。Roary的输入格式就是prokka生成的那种GFF3格式它会从GFF文件里提取每个样本的基因信息然后根据序列相似性做直系同源聚类最终生成泛基因组分析结果。# Roary对prokka结果的直接消费无需额外转换 roary -f roary_output -e -n -v *.gff我自己做过多株大肠杆菌的泛基因组分析流程就是prokka → Roary → 下游可视化。如果不用prokka而用其他工具后续还需要写脚本转换格式非常麻烦。6.2 基因簇注释与次级代谢产物预测做链霉菌、芽孢杆菌这类次级代谢产物丰富的菌株时我通常会在prokka注释完成之后直接把.gbk文件丢给antiSMASH预测次级代谢基因簇。antiSMASH对输入的格式要求是GenBank或EMBL格式prokka生成的.gbk和.sqn都可以直接作为输入。不过需要注意一点prokka注释的是整个基因组的基因结构antiSMASH做的是基于基因簇保守模块的识别两者在边界判断上偶尔会有出入。遇到antiSMASH预测出一个非常大的基因簇而prokka只注释出中间几个CDS时通常是prokka漏掉了边界基因。这时候需要手动检查该区域的序列补充注释后再进下游分析。6.3 系统发育分析中的应用prokka的另一个高频应用场景是提供系统发育分析用的直系同源蛋白。将多个亲缘关系较近菌株的prokka蛋白序列通过OrthoFinder或者Roary聚类后提取单拷贝直系同源基因串联比对后建树。这个方法的优势在于prokka统一的标准输出保证了不同样本之间基因命名的可比性避免了不同注释工具产出的差异变成伪信号。6.4 关于“真核项目”数据的补充说明前面提到prokka不适合真核基因组注释但实际工作中总会遇到手里拿着真核数据的情况。比如最近不少人从公共数据库下载植物基因组包括芍药这类园艺植物的基因组文件和注释文件。这类文件通常是FASTA格式的基因组序列加GFF3格式的注释文件下载后应该先验证注释文件与基因组序列的染色体命名是否一致再决定后续分析。不少人在这一步翻车就是因为GFF里的染色体ID和FASTA里的头部名称对不上下游所有分析全部错乱。这个方法论虽然在植物基因组场景下更常见但对于养成“先验证再分析”的习惯非常有帮助——真核也好、原核也好多一步验证后面少十步返工。结尾的话我在实际项目里大规模用过prokka之后才体会到真正好的生信工具不是功能越多越好而是把用户从无谓的格式转换和脚本编写中解放出来让你把精力放到生物学问题上。prokka就是这样的工具——它未必每个环节都做到最好但整个流程用得顺手结果可信和上下游工具的衔接做得漂亮。最后分享一个小习惯给大家我每次跑完prokka从来不会直接删掉annotation目录。除了常规的.gff和.gbk之外subtilis.log文件值得永久保留。这个日志文件里记录了所用的prokka版本号、参数组合、运行时间和输入文件的基本统计信息。等到下个月你回想起来“咦上次那个注释是怎么跑的”的时候这个日志文件就是最快的线索比任何笔记都靠谱。