做二代全基因组测序的样本里线粒体基因组几乎是最容易顺手拿到的副产品也是最好发的数据。但真到自己上手跑的时候很多人会卡在一个很尴尬的位置数据在硬盘里躺着MitoZ 也装了命令抄下来了跑完却不知道结果到底能不能用。我前后在鱼类、昆虫、贝类几类样本上跑过上百次踩过的坑基本覆盖了这个工具从安装到出图的所有环节。这篇就把动物线粒体基因组组装注释这条流水线拆开讲透重点不在复述说明书而在那些说明书上不会写、但你一定会撞上的问题。MitoZ 是一个把线粒体基因组组装、注释、可视化打包成一条命令的工具输入是双端短读长原始数据输出是一条或几条线粒体序列、一份 GenBank 注释文件和一张环形图。它适合三类人手里有 WGS 数据想顺带把线粒体拿出来的人、做系统发育需要大量线粒体序列的人、以及刚接触组装流程想先跑通一个完整案例的新手。它对使用者的要求不高会敲命令行就行但它对输入数据质量和参数匹配度的要求很高这是它最大的特点也是后面所有坑的来源。1. MitoZ 到底替我们做了哪几件事1.1 五个环节串成的一条流水线把 MitoZ 当成黑箱用是最容易出问题的用法。它的内部其实是五个相对独立的环节每个环节都有自己的失败模式知道它卡在哪一步排查效率会差好几倍。第一步读长过滤。调用 fastp 对原始双端数据做质控去掉接头、低质量碱基和过短的读长。这一步本身没什么技术含量但它决定了后面组装的上限。第二步从头组装。默认调 MEGAHIT也可以切到另一套组装器。注意这里组的是全部数据不是只组线粒体——MitoZ 没有本事在组装前把线粒体 reads 挑出来它只能在组装后筛。第三步钓取候选序列。把组装出来的所有 contig 拿去和内置的线粒体基因数据库做比对蛋白编码基因用翻译比对rRNA 用核酸比对有命中的 contig 被标记为候选。第四步回收重组装。这是 MitoZ 相比组完直接 BLAST 挑更有价值的地方把比对到候选 contig 的 reads 重新捞出来再做一轮组装绕开核基因组 reads 的干扰把线粒体序列拼得更完整。第五步注释与出图。tRNA 用 tRNAscan-SE 扫rRNA 和蛋白编码基因用比对加结构域搜索确定边界最后调 Circos 画环图、做基因顺序比较。这条链路里第二、三、四步是决定成败的核心第五步只影响结果好不好看、准不准。新手最常见的误区是把注意力全放在最后那张图上其实图上缺个 tRNA问题往往出在第四步拼出来的序列本来就是断的。1.2 一次完整的计算到底花在哪很多人第一次跑会惊讶于耗时——明明线粒体只有十几 kb为什么跑了好几个小时。原因很简单组装那一步处理的是全基因组数据。以 5 Gb 的双端 150 bp 数据为例大约 1670 万对 reads。MEGAHIT 在这种量级上普通服务器 20 线程、64 GB 内存组装本身就要一到两个小时后面的比对和回收再叠加半小时到一小时。真正只跟线粒体有关的计算可能只占总时间的 10%。理解了这一点很多优化思路就自然出来了如果你只关心线粒体完全没必要用全部数据。我通常的做法是先做一次数据抽样把 reads 随机抽到 12 Gb 再喂给 MitoZ跑通流程、看看结果质量确认没问题后再决定要不要用全量数据重跑。线粒体在总 DNA 里的占比通常在 0.5%10% 之间肌肉、肝脏这类组织偏高血液、毛发、陈旧标本偏低抽到 1 Gb 时线粒体 reads 仍然有几十万条覆盖度早就过百了。这里可以算一笔账。假设抽 1 Gb 数据线粒体占比保守按 0.5% 算来自线粒体的碱基数就是 5 Mb除以 16.5 kb 的线粒体基因组长度覆盖度约 300×。这个深度对组装来说绰绰有余甚至偏冗余。真正的瓶颈从来不是深度不够而是核基因组里的线粒体假基因numt在捣乱这个后面单独讲。1.3 和 GetOrganelle、NOVOPlasty 这类工具怎么选同类工具不止一个选错了会浪费大量时间。我把几个常见方案的差异整理成一张表判断依据就三条有没有近缘参考、数据量多少、要不要一步出注释。工具是否需要种子序列数据量要求注释能力适用场景MitoZ不需要自带数据库自动钓取中等110 Gb内置注释一条命令出 GenBank类群常见、想一步到位GetOrganelle需要近缘参考或内置库较宽松需要另配注释工具有可靠近缘参考、追求组装质量NOVOPlasty需要一个种子序列如 cox1 片段较小无数据量紧张、有明确种子纯手工组装比对手工注释需要参考任意全手工类群特殊、自动注释老出错我的实际选择逻辑是先看类群。脊椎动物、节肢动物、软体动物这些线粒体研究充分的类群MitoZ 的命中率很高直接上就行。如果是某些冷门无脊椎动物内置数据库覆盖不全MitoZ 可能钓不到东西这时候换 GetOrganelle 或者干脆手工从总基因组组装结果里把线粒体 contig 抠出来更靠谱。再看数据形态。MitoZ 只吃二代短读长。如果手上是三代长读长数据MitoZ 帮不上忙得先用长读长专用组装器拼出环状序列再单独调annotate模块做注释。这一点很多人搞混把三代数据硬塞进去结果当然是失败的。最后提醒一句这几个工具的结果可以互相验证但不能互相替代。我经常在 MitoZ 结果不理想时用 GetOrganelle 再跑一遍两条序列比对一下差异区域往往就是真正的难点所在。2. 环境搭建90% 的失败都发生在这一步2.1 依赖清单与版本地雷MitoZ 不是一个轻量脚本它是一堆第三方工具的上层封装。一份完整的运行环境大致需要这些组件Python 运行环境不同大版本对 Python 的要求不一样这一点必须先确认装错了后面全白搭组装器MEGAHIT或另一套基于 de Bruijn 图的组装器序列比对BLAST 工具集包含核酸比对和翻译比对两个程序结构域搜索HMMER 套件tRNA 识别tRNAscan-SE多序列比对用于基因边界确认和序列对齐绘图Circos依赖 Perl 及其模块、R 及其绘图包质控fastp版本地雷主要有三处我按踩坑频率排序第一处BLAST 版本。新版 BLAST 在输出格式和参数默认值上有调整老脚本里的一些写法会失效表现为比对结果为空或者程序直接报错。如果你发现流程卡在比对环节、日志里什么都没留下先查 BLAST 版本。第二处Python 版本。早期版本依赖 Python 2.7新版迁到了 Python 3。混用会报语法错误那种SyntaxError出现在第三方库文件里的报错基本都是这个原因。装之前先看版本对应的 Python 需求不要想当然。第三处Perl 与 Circos。这个只在最后出图时暴露出来表现为前面几个小时全跑完了最后一步莫名其妙挂掉。Circos 的 Perl 模块依赖非常琐碎缺一个就报错而且报错信息指向性很差。提示如果你的样本量不大优先考虑用官方打包好的完整分发包或者容器镜像跳过所有依赖编译这一步能省掉大半天。2.2 三条安装路线的取舍实际可选的就三条路各有适用场景。路线一Conda 环境安装。最省事一行命令拉下来。优点是干净、可复现、删了不留痕迹。缺点是如果某个依赖在软件源里的版本没对齐跑起来会在某个隐蔽的地方崩掉排查起来很痛苦。适合已经熟悉这套工具链的人。路线二官方打包分发包。作者把所有依赖预编译好打成一个压缩包解压即用。优点是版本组合是经过测试的稳定性最高缺点是对系统底层库有要求跨发行版可能跑不起来而且路径写死了不好挪。适合第一次上手、想先看到结果的人。路线三容器镜像。隔离性最好跟宿主机环境完全解耦。缺点是数据目录挂载、线程数、内存限制这些参数要额外配跑大数据时容易因为容器资源限制而失败。适合在共享服务器上跑、不想污染环境的情况。我自己的习惯是先容器跑通、验证流程可用再在 Conda 里搭一套长期用的。这样既保证了第一次能出结果又保证后面批量跑的时候效率高。2.3 装完必须跑的六条自检命令装完直接上真实数据是典型的错误做法。我给自己定的规矩是先跑六条检查全部通过再动数据主程序能否正常打印帮助信息——验证入口和 Python 环境组装器能否打印版本——验证组装环节核酸比对程序能否打印版本——验证比对环节翻译比对程序能否打印版本——这一步单独查因为两个程序有时来自不同安装源结构域搜索程序能否打印版本——验证注释环节tRNA 识别程序能否打印版本——验证 tRNA 注释第 4 条最容易被跳过。我遇到过一次核酸比对能跑翻译比对找不到结果就是 tRNA 和 rRNA 都注释出来了13 个蛋白编码基因一个都没有。日志里只有一行比对失败没有任何提示。还有一个隐藏点如果某个依赖不在系统 PATH 里大部分模块支持用参数手动指定路径比如 fastp。多环境共存的时候这个参数能救命值得提前记下来。3. 参数不是玄学每个开关背后的判断依据3.1 clade 与 geneticcode 必须成对确定这两个参数是整条流程里最容易出错、后果最严重的两个。它们决定了用哪套遗传密码去翻译序列、用哪套数据库去钓取基因。clade告诉程序你的物种属于哪个大类群从而选择对应的参考数据库。填错了钓取阶段就找不到东西后面的组装回收也无从谈起。geneticcode告诉程序用哪张密码子表翻译。填错了序列拼出来了但蛋白编码基因会被翻译成一堆提前终止的乱码注释结果惨不忍睹。常见的对应关系我整理如下具体到目、科一级最好还是查一下该类的文献惯例类群建议 clade 值建议 geneticcode脊椎动物鱼、两栖、爬行、鸟、哺乳Chordata2节肢动物昆虫、甲壳、蛛形Arthropoda5软体动物贝、螺、头足Mollusca5环节动物Annelida5线虫Nematoda5棘皮动物Echinodermata9扁形动物Platyhelminthes9刺胞动物Cnidaria多为 1 或 4按类群核对判断有没有填错有个很直接的信号看注释出来的蛋白编码基因是不是普遍偏短、且每条都带提前终止密码子。如果是八成是密码子表错了。另一种信号是 13 个基因只注释出三四个那可能是 clade 错了导致钓取失败。注意clade 的取值在不同版本里可能有增减敲帮助命令看当前版本支持哪些值别照抄网上的老帖子。3.2 assembler、内存与数据量的三角关系组装器的选择不是随便挑一个它跟你的数据量和机器内存强相关。默认的组装器MEGAHIT内存效率高、速度快是绝大多数情况下的首选。它的设计目标就是用有限内存处理海量数据几 Gb 的数据在几十 GB 内存的机器上能跑得很舒服。我一般优先用它。另一套组装器在某些高杂合或者特殊结构的样本上表现可能更好但代价是内存占用高得多且速度慢。只有在默认组装器拼出来的结果明显有问题比如线粒体被拆成很多短片段时才值得切过去试。内存是硬约束不是建议值。我见过最典型的失败场景是跑到组装阶段进程被系统直接杀掉日志里什么都没有只在系统日志里留下一行内存不足。这种情况你调多少线程都没用只能降数据量或者加内存。我常用的做法是把数据抽到组装器舒适区如果机器只有 32 GB 内存就把输入数据控制在 2 Gb 以内有 128 GB 就把上限提到 810 Gb。线粒体覆盖度早就够了多喂的数据只是在浪费时间。3.3 run_mode、dbnum、thread_number 的取舍这三个参数属于调了有感觉但不致命的类型理解它们能帮你优化时间不理解也不影响正确性。run_mode控制流程的彻底程度。快速模式会跳过一些耗时的搜索步骤适合先摸底、看这批数据有没有希望完整模式会做多轮钓取和回收结果更完整但慢。我的一般策略是先用快速模式跑一批样本看整体成功率对快速模式结果不理想的样本再单独用完整模式重跑。这样避免了所有样本都跑完整流程的时间浪费。dbnum控制把参考数据库切成几份并行搜索。这个参数跟线程数配合使用——数据库切得太少线程跑不满切得太多每个分片太小比对本身的固定开销占比上升。经验值是切成线程数的四分之一到二分之一比如 20 线程切成 510 份。thread_number就是线程数没什么好说的但有个细节这个值要同时满足组装器和比对程序的需求有些环境里它不能超过物理核心数超了反而变慢。3.4 两条可以直接抄走的命令说完参数给两条我自己常用的命令模板。第一条是完整的从头跑mitoz All \ --fq1 sample_R1.fq.gz \ --fq2 sample_R2.fq.gz \ --outprefix sample_mito \ --clade Chordata \ --geneticcode 2 \ --thread_number 20 \ --dbnum 8 \ --max_read_length 150 \ --run_mode 2 \ --fastp /opt/biosoft/fastp第二条是已经有序列、只想注释的情况——这是很实用的场景比如你自己手工拼好了环状序列或者从总基因组组装结果里抠出了一条 contigmitoz annotate \ --genefasta candidate_mito.fa \ --outprefix sample_mito \ --clade Arthropoda \ --geneticcode 5 \ --thread_number 10第二条命令省掉了组装环节几分钟就能出结果非常适合做参数试错。不确定 clade 或 geneticcode 填什么的时候先用第二条命令拿同一条序列反复试看哪组参数注释出的基因最全、最完整再去跑完整流程。这个技巧帮我节省了大量时间。4. 结果怎么读哪些文件有用哪些是噪音4.1 输出目录速查表一次完整运行会在输出前缀目录下生成一大堆文件第一次看会懵。实际需要关注的只有几个文件/目录内容重要程度前缀.fa最终组装的线粒体序列最关键前缀.gbGenBank 格式注释文件最关键前缀.summary.txt序列与注释的汇总统计高前缀.circos.pdf / .svg环形图中出图用前缀.gene.order.png基因顺序比较图中验证用前缀.depth.stat覆盖度统计高判断质量用前缀.trna.txt / .rrna.txt / .gene.txt各类基因的分项结果高核对用中间目录质控、组装、比对、注释、绘图各阶段中间文件低排查时才看我的习惯是先看.fa的长度再看.gb里基因数然后看覆盖度统计最后才看图。这个顺序能最快判断这条序列值不值得继续投入时间。中间目录不要急着删。排查问题时它们就是现场证据尤其是组装目录和比对目录能告诉你钓取这一步到底有没有命中以及命中了什么。4.2 判断组装质量的四个硬指标拿到.fa之后用四个指标快速打分基本能判断这条序列能不能用。指标一总长度是否落在该类群的正常区间。脊椎动物多在 16.5 kb 上下昆虫多在 1517 kb线虫偏短1314 kb软体动物变异最大从 14 kb 到 20 kb 都有。如果拼出来只有 8 kb那肯定是断了如果拼出 25 kb那可能混进了重复序列或者 numt。指标二覆盖度是否均匀。打开.depth.stat看沿序列的覆盖度分布。一条健康的线粒体序列覆盖度应该是平缓的高原偶尔在控制区有个小起伏控制区有串联重复比对时容易堆叠。如果出现一段覆盖度突然掉到个位数甚至 0说明那里有断点或者组装错误。指标三是否能成环。真正的线粒体基因组是环状的。判断方法是看序列两端是否存在重叠——如果末端几十到几百个碱基和开头高度相似那基本就是环状结构被线性化时的产物手动切掉重叠加起来就行。反过来如果两端完全对不上、覆盖度也不衔接那就是断的需要补洞。指标四基因数是否完整。动物线粒体基因组成非常保守13 个蛋白编码基因 22 个 tRNA 2 个 rRNA 37 个。少一两个可能是注释灵敏度问题比如某个 tRNA 反密码子特殊少五六个就是序列本身的问题了。4.3 注释结果的核对逻辑13 22 2数量对上了不代表质量对上了还得看细节。我核对的顺序是固定的先看蛋白编码基因。逐个检查三件事长度是否合理、起始密码子是否正常、有没有提前终止。有一个例外必须记住——有些基因在序列上确实会出现提前终止或者不完整末端因为它们的转录本要经过多聚腺苷化来补齐这是线粒体的正常现象不是错误。判断方法是把序列翻译出来看是不是只有一个位置异常、其余部分完全连贯。如果断得零零碎碎那就是真错了。再看 tRNA。22 个 tRNA 覆盖 20 种氨基酸其中亮氨酸和丝氨酸各有两个反密码子不同。少了某个 tRNA先别急着改序列可以先换一个 tRNA 识别工具交叉验证。线粒体 tRNA 二级结构短、变异大单一工具漏检很常见。低灵敏度模式或者换个算法往往就能捡回来。最后看 rRNA。两个 rRNA大亚基和小亚基的边界靠比对和二级结构确定自动注释经常有几十个碱基的偏差。这个偏差在系统发育分析里影响不大但如果要提交数据库最好手工修一下。还有一个必查项基因之间的重叠与间隔。动物线粒体里有几个固定的重叠模式比如相邻的两个 ATP 合酶亚基基因常重叠几个碱基、两个 NADH 脱氢酶亚基基因也常重叠。这些重叠如果在你的注释结果里消失了说明边界定错了。反过来如果出现大片不该有的基因间区也可能是边界问题或者拼装错误。5. 避坑实录报错、烂结果与补救办法5.1 运行期报错的排查顺序报错不可怕可怕的是没日志。我遇到过的报错按频率排下来大致是这几类排查顺序也有讲究。第一类找不到命令或者依赖。报错信息里会出现command not found之类的字样。这是环境问题八成是依赖没装全或者不在 PATH 里。解决办法是用参数显式指定路径或者把依赖目录加进环境变量。第二类内存不足被系统杀掉。特征是没有报错、进程莫名消失、退出码异常。查系统日志能看到内存耗尽的记录。对策是降数据量或者加内存别指望调参数能绕过去。第三类比对或搜索返回空结果。这类最麻烦因为流程可能正常结束但什么都没注释出来。原因可能是 clade 填错、参考库版本不对、序列本身不是线粒体。排查方法是去中间目录里翻比对输出文件看是不是真的一条命中都没有。第四类绘图阶段崩溃。前面全跑完最后一步挂掉。九成是绘图软件的环境问题Perl 模块缺失、R 包缺失。这类问题不影响数据本身.fa和.gb都是完整的可以手工用别的工具出图。心得养成先看中间目录再问人的习惯。一个跑失败的流程中间目录里留下的信息量远大于一句报错了。5.2 组装出来一堆 contig 怎么办这是最典型的结果不佳场景.fa里躺着七八条 contig长短不一加起来才勉强够线粒体长度。原因就那么几个我按可能性排序。原因一核基因组里的线粒体假基因numt干扰。这是头号杀手。动物的核基因组里普遍插入了线粒体片段的拷贝这些片段和真线粒体高度相似组装时会互相抢reads导致线粒体被拆成多个片段甚至拼出嵌合序列。怎么识别 numt三个特征覆盖度明显偏低numt 的覆盖度等于核基因组覆盖度通常比线粒体低一到两个数量级、含有提前终止密码子或移码因为插入后不再受纯化选择约束会积累突变、碱基组成偏离真线粒体通常 AT 含量很高numt 会向核基因组的 GC 水平漂移。用 reads 回贴到每条 contig 上看覆盖度一眼就能分出来。原因二数据量不够。特别是保存久的标本、非损伤取样如粪便、毛发的样本线粒体 reads 本来就少。这种情况下补测序是唯一解靠调参数没用。原因三控制区太难拼。控制区富含串联重复长度变异大是所有组装器的共同噩梦。这里我要说一句可能反直觉的话控制区拼不完整很多时候不影响使用。如果做的是系统发育或者基因组成分析把控制区标注为不完整区域即可不必死磕。原因四样本里混了别的物种。比如寄生虫、共生菌、食物残留。这类序列的 GC 含量和密码子偏好跟目标物种差得很远比对一下就能发现。对策是在组装前做去宿主或者按 GC/覆盖度筛 contig。补救的做法把目标 contig 挑出来用 reads 回贴看覆盖度把明显是 numt 的剔除然后对剩下的 contig 做延长拼接或者设计引物做 PCR 补洞。MitoZ 的注释模块可以单独对拼好的序列再跑一遍不用从头重来。5.3 注释缺基因、多基因、边界错位缺基因最常见的是缺 tRNA 和缺那个最短的 ATP 合酶亚基基因。前者是识别灵敏度问题换工具或者降低阈值能救后者是因为它实在太短很多类群只有五六十个氨基酸比对打分上不去被过滤掉了。对策是手工按位置和保守区找或者降低比对阈值重跑。多基因的情况通常是同一个基因被注释了两遍——典型原因是序列里存在重复片段numt 残留或者组装错误导致某段重复或者基因重叠区域被重复计算。看到 38 个、39 个基因先查重复。边界错位最隐蔽。表现为基因长度多了或少了十几到几十个碱基起始密码子不对。这种情况一般不影响下游分析但如果要提交数据库还是手工修一下比较稳妥。修的依据是翻译成蛋白后和近缘物种的同源蛋白比对看哪个位置开始出现连续的错配。还有一个特殊情况值得单独说基因重排。某些类群的线粒体基因顺序和标准顺序差别很大比如部分软体动物、部分寄生虫。如果注释结果里基因顺序乱得离谱先别怀疑工具去查一下这个类群是不是本来就有重排。这类信息在基因顺序比较图上能看得很清楚。5.4 常见问题速查表把上面这些浓缩成一张表出问题时按症状查症状最可能的原因处理办法流程跑完但注释出 0 个基因clade 填错或参考库不匹配换 clade 用注释模块单独重试蛋白编码基因普遍带提前终止geneticcode 填错按类群换密码子表重跑注释组装结果为多条短 contignumt 干扰或数据量不足按覆盖度筛 contig必要时补测序某段覆盖度骤降到接近 0组装断点用 reads 回贴定位设计引物补洞总长度远超正常范围混入重复序列或嵌合检查末端重叠切掉冗余部分tRNA 数量少于 22识别灵敏度不足换工具交叉验证降低阈值缺少最短的 ATP 合酶基因序列太短打分不够按保守区手工定位绘图阶段崩溃绘图环境缺失单独用其他工具出图不影响数据跑完进程莫名消失内存不足降数据量或加内存6. 收尾从自动结果到能提交的序列6.1 断点定位与补洞的实际做法自动流程给出的序列只能算半成品。要变成一条干净、完整、能提交的序列还需要几步手工活。第一步确认断点位置。把原始 reads 回贴到组装序列上看覆盖度在哪里断掉。断点通常表现为覆盖度急剧下降或者出现一堆 soft-clip 的 reads。把断点前后各一两百个碱基记下来这就是需要补的区域。第二步尝试延长。很多时候断点并不是真的没 data而是组装器在那里放弃了。可以把断点两端的序列当锚点去原始 reads 里捞跨越断点的 reads做局部延长。或者用另一套组装器重跑一遍两套结果互相补。第三步PCR 补洞。如果数据里确实没有跨越断点的 reads通常是重复区或者极端 AT 富集区那就只能回到实验台。在断点两侧设计一对引物扩增出产物做 Sanger 测序把序列接上去。这一步听起来麻烦但一条 16 kb 的线粒体通常也就一两个洞一两天能搞定。第四步确定起点和方向。数据库提交和后续分析通常要求统一的起点。常见约定是以某个特定基因或者 tRNA 作为序列开头。取哪个起点看你的研究类群惯例也可以参考近缘物种已提交的序列。定好之后把环状序列旋转一下切掉末端重叠得到一条标准的线性表示。提示补洞前的序列和补洞后的序列一定要都留着比对一下确认改动区域没有引入新的错误。我有一次补完洞忘了核对方向结果整条序列反向互补了白折腾半天。6.2 提交前的自查清单序列收拾干净之后提交之前我会过一遍这几项每一项都对应过真实的返工长度是否在类群正常区间是否确认为环状末端重叠已处理基因数是否为 13 22 2每个蛋白编码基因的起始、终止密码子是否合理是否存在异常的大段基因间区超过几百个碱基的空白需要解释序列的碱基组成是否符合该类群特征AT 含量、GC skew控制区是否标注清楚不完整是否说明注释文件的格式是否符合数据库要求是否有 numt 嫌疑片段残留其中异常大段基因间区这一项最容易被忽略。动物线粒体基因排布非常紧凑基因间区通常只有几十个碱基。如果你发现某两个基因之间空了两千个碱基那几乎可以肯定是要么拼错了要么混进了外源序列要么那段是没注释出来的基因。6.3 出图与二次美化自动生成的环形图够用但正式发表或者做汇报时往往需要再加工。我的经验是分工明确用自动工具出准确的图用绘图工具出好看的图。如果自动环形图能正常出直接拿矢量格式PDF 或 SVG进矢量编辑软件微调改字体、调颜色、补图例比重新画一遍快得多。如果自动出图失败可以退一步用在线注释工具重跑注释再出图或者自己用绘图代码画——把注释结果里每个基因的起止坐标读出来画成同心圆弧工作量不大且完全可控。基因顺序比较图的价值容易被低估。它不只是一张配图在判断组装和注释是否正确时非常有用把目标物种和几个近缘物种的基因顺序排在一起如果目标物种出现了大段倒位或者易位要么是真的有重排值得写进文章要么是拼装或注释出了问题需要回头查。两种情况都需要你重视。最后分享一个我自己常用的判断习惯拿到任何一条自动组装的线粒体序列先把 13 个蛋白编码基因翻译出来把氨基酸序列拼在一起看一遍。如果中间出现了莫名其妙的杂合、连续错配、或者某一段明显不像同源序列那问题就在那段。这个方法比看任何统计指标都直观几分钟就能定位到具体位置。另外工具本身在持续更新参数名和默认行为会变。我现在的做法是把每次跑通的命令和对应的版本号记在一个文本文件里下次直接复制。这个习惯帮我省掉过无数次上次明明能跑的困惑——因为上次跑的确实不是同一个版本。 
长河君 · 2026/11/3 5:24:16 · 阅读 1.2k MitoZ动物线粒体基因组组装注释全流程与避坑指南
做二代全基因组测序的样本里线粒体基因组几乎是最容易顺手拿到的副产品也是最好发的数据。但真到自己上手跑的时候很多人会卡在一个很尴尬的位置数据在硬盘里躺着MitoZ 也装了命令抄下来了跑完却不知道结果到底能不能用。我前后在鱼类、昆虫、贝类几类样本上跑过上百次踩过的坑基本覆盖了这个工具从安装到出图的所有环节。这篇就把动物线粒体基因组组装注释这条流水线拆开讲透重点不在复述说明书而在那些说明书上不会写、但你一定会撞上的问题。MitoZ 是一个把线粒体基因组组装、注释、可视化打包成一条命令的工具输入是双端短读长原始数据输出是一条或几条线粒体序列、一份 GenBank 注释文件和一张环形图。它适合三类人手里有 WGS 数据想顺带把线粒体拿出来的人、做系统发育需要大量线粒体序列的人、以及刚接触组装流程想先跑通一个完整案例的新手。它对使用者的要求不高会敲命令行就行但它对输入数据质量和参数匹配度的要求很高这是它最大的特点也是后面所有坑的来源。1. MitoZ 到底替我们做了哪几件事1.1 五个环节串成的一条流水线把 MitoZ 当成黑箱用是最容易出问题的用法。它的内部其实是五个相对独立的环节每个环节都有自己的失败模式知道它卡在哪一步排查效率会差好几倍。第一步读长过滤。调用 fastp 对原始双端数据做质控去掉接头、低质量碱基和过短的读长。这一步本身没什么技术含量但它决定了后面组装的上限。第二步从头组装。默认调 MEGAHIT也可以切到另一套组装器。注意这里组的是全部数据不是只组线粒体——MitoZ 没有本事在组装前把线粒体 reads 挑出来它只能在组装后筛。第三步钓取候选序列。把组装出来的所有 contig 拿去和内置的线粒体基因数据库做比对蛋白编码基因用翻译比对rRNA 用核酸比对有命中的 contig 被标记为候选。第四步回收重组装。这是 MitoZ 相比组完直接 BLAST 挑更有价值的地方把比对到候选 contig 的 reads 重新捞出来再做一轮组装绕开核基因组 reads 的干扰把线粒体序列拼得更完整。第五步注释与出图。tRNA 用 tRNAscan-SE 扫rRNA 和蛋白编码基因用比对加结构域搜索确定边界最后调 Circos 画环图、做基因顺序比较。这条链路里第二、三、四步是决定成败的核心第五步只影响结果好不好看、准不准。新手最常见的误区是把注意力全放在最后那张图上其实图上缺个 tRNA问题往往出在第四步拼出来的序列本来就是断的。1.2 一次完整的计算到底花在哪很多人第一次跑会惊讶于耗时——明明线粒体只有十几 kb为什么跑了好几个小时。原因很简单组装那一步处理的是全基因组数据。以 5 Gb 的双端 150 bp 数据为例大约 1670 万对 reads。MEGAHIT 在这种量级上普通服务器 20 线程、64 GB 内存组装本身就要一到两个小时后面的比对和回收再叠加半小时到一小时。真正只跟线粒体有关的计算可能只占总时间的 10%。理解了这一点很多优化思路就自然出来了如果你只关心线粒体完全没必要用全部数据。我通常的做法是先做一次数据抽样把 reads 随机抽到 12 Gb 再喂给 MitoZ跑通流程、看看结果质量确认没问题后再决定要不要用全量数据重跑。线粒体在总 DNA 里的占比通常在 0.5%10% 之间肌肉、肝脏这类组织偏高血液、毛发、陈旧标本偏低抽到 1 Gb 时线粒体 reads 仍然有几十万条覆盖度早就过百了。这里可以算一笔账。假设抽 1 Gb 数据线粒体占比保守按 0.5% 算来自线粒体的碱基数就是 5 Mb除以 16.5 kb 的线粒体基因组长度覆盖度约 300×。这个深度对组装来说绰绰有余甚至偏冗余。真正的瓶颈从来不是深度不够而是核基因组里的线粒体假基因numt在捣乱这个后面单独讲。1.3 和 GetOrganelle、NOVOPlasty 这类工具怎么选同类工具不止一个选错了会浪费大量时间。我把几个常见方案的差异整理成一张表判断依据就三条有没有近缘参考、数据量多少、要不要一步出注释。工具是否需要种子序列数据量要求注释能力适用场景MitoZ不需要自带数据库自动钓取中等110 Gb内置注释一条命令出 GenBank类群常见、想一步到位GetOrganelle需要近缘参考或内置库较宽松需要另配注释工具有可靠近缘参考、追求组装质量NOVOPlasty需要一个种子序列如 cox1 片段较小无数据量紧张、有明确种子纯手工组装比对手工注释需要参考任意全手工类群特殊、自动注释老出错我的实际选择逻辑是先看类群。脊椎动物、节肢动物、软体动物这些线粒体研究充分的类群MitoZ 的命中率很高直接上就行。如果是某些冷门无脊椎动物内置数据库覆盖不全MitoZ 可能钓不到东西这时候换 GetOrganelle 或者干脆手工从总基因组组装结果里把线粒体 contig 抠出来更靠谱。再看数据形态。MitoZ 只吃二代短读长。如果手上是三代长读长数据MitoZ 帮不上忙得先用长读长专用组装器拼出环状序列再单独调annotate模块做注释。这一点很多人搞混把三代数据硬塞进去结果当然是失败的。最后提醒一句这几个工具的结果可以互相验证但不能互相替代。我经常在 MitoZ 结果不理想时用 GetOrganelle 再跑一遍两条序列比对一下差异区域往往就是真正的难点所在。2. 环境搭建90% 的失败都发生在这一步2.1 依赖清单与版本地雷MitoZ 不是一个轻量脚本它是一堆第三方工具的上层封装。一份完整的运行环境大致需要这些组件Python 运行环境不同大版本对 Python 的要求不一样这一点必须先确认装错了后面全白搭组装器MEGAHIT或另一套基于 de Bruijn 图的组装器序列比对BLAST 工具集包含核酸比对和翻译比对两个程序结构域搜索HMMER 套件tRNA 识别tRNAscan-SE多序列比对用于基因边界确认和序列对齐绘图Circos依赖 Perl 及其模块、R 及其绘图包质控fastp版本地雷主要有三处我按踩坑频率排序第一处BLAST 版本。新版 BLAST 在输出格式和参数默认值上有调整老脚本里的一些写法会失效表现为比对结果为空或者程序直接报错。如果你发现流程卡在比对环节、日志里什么都没留下先查 BLAST 版本。第二处Python 版本。早期版本依赖 Python 2.7新版迁到了 Python 3。混用会报语法错误那种SyntaxError出现在第三方库文件里的报错基本都是这个原因。装之前先看版本对应的 Python 需求不要想当然。第三处Perl 与 Circos。这个只在最后出图时暴露出来表现为前面几个小时全跑完了最后一步莫名其妙挂掉。Circos 的 Perl 模块依赖非常琐碎缺一个就报错而且报错信息指向性很差。提示如果你的样本量不大优先考虑用官方打包好的完整分发包或者容器镜像跳过所有依赖编译这一步能省掉大半天。2.2 三条安装路线的取舍实际可选的就三条路各有适用场景。路线一Conda 环境安装。最省事一行命令拉下来。优点是干净、可复现、删了不留痕迹。缺点是如果某个依赖在软件源里的版本没对齐跑起来会在某个隐蔽的地方崩掉排查起来很痛苦。适合已经熟悉这套工具链的人。路线二官方打包分发包。作者把所有依赖预编译好打成一个压缩包解压即用。优点是版本组合是经过测试的稳定性最高缺点是对系统底层库有要求跨发行版可能跑不起来而且路径写死了不好挪。适合第一次上手、想先看到结果的人。路线三容器镜像。隔离性最好跟宿主机环境完全解耦。缺点是数据目录挂载、线程数、内存限制这些参数要额外配跑大数据时容易因为容器资源限制而失败。适合在共享服务器上跑、不想污染环境的情况。我自己的习惯是先容器跑通、验证流程可用再在 Conda 里搭一套长期用的。这样既保证了第一次能出结果又保证后面批量跑的时候效率高。2.3 装完必须跑的六条自检命令装完直接上真实数据是典型的错误做法。我给自己定的规矩是先跑六条检查全部通过再动数据主程序能否正常打印帮助信息——验证入口和 Python 环境组装器能否打印版本——验证组装环节核酸比对程序能否打印版本——验证比对环节翻译比对程序能否打印版本——这一步单独查因为两个程序有时来自不同安装源结构域搜索程序能否打印版本——验证注释环节tRNA 识别程序能否打印版本——验证 tRNA 注释第 4 条最容易被跳过。我遇到过一次核酸比对能跑翻译比对找不到结果就是 tRNA 和 rRNA 都注释出来了13 个蛋白编码基因一个都没有。日志里只有一行比对失败没有任何提示。还有一个隐藏点如果某个依赖不在系统 PATH 里大部分模块支持用参数手动指定路径比如 fastp。多环境共存的时候这个参数能救命值得提前记下来。3. 参数不是玄学每个开关背后的判断依据3.1 clade 与 geneticcode 必须成对确定这两个参数是整条流程里最容易出错、后果最严重的两个。它们决定了用哪套遗传密码去翻译序列、用哪套数据库去钓取基因。clade告诉程序你的物种属于哪个大类群从而选择对应的参考数据库。填错了钓取阶段就找不到东西后面的组装回收也无从谈起。geneticcode告诉程序用哪张密码子表翻译。填错了序列拼出来了但蛋白编码基因会被翻译成一堆提前终止的乱码注释结果惨不忍睹。常见的对应关系我整理如下具体到目、科一级最好还是查一下该类的文献惯例类群建议 clade 值建议 geneticcode脊椎动物鱼、两栖、爬行、鸟、哺乳Chordata2节肢动物昆虫、甲壳、蛛形Arthropoda5软体动物贝、螺、头足Mollusca5环节动物Annelida5线虫Nematoda5棘皮动物Echinodermata9扁形动物Platyhelminthes9刺胞动物Cnidaria多为 1 或 4按类群核对判断有没有填错有个很直接的信号看注释出来的蛋白编码基因是不是普遍偏短、且每条都带提前终止密码子。如果是八成是密码子表错了。另一种信号是 13 个基因只注释出三四个那可能是 clade 错了导致钓取失败。注意clade 的取值在不同版本里可能有增减敲帮助命令看当前版本支持哪些值别照抄网上的老帖子。3.2 assembler、内存与数据量的三角关系组装器的选择不是随便挑一个它跟你的数据量和机器内存强相关。默认的组装器MEGAHIT内存效率高、速度快是绝大多数情况下的首选。它的设计目标就是用有限内存处理海量数据几 Gb 的数据在几十 GB 内存的机器上能跑得很舒服。我一般优先用它。另一套组装器在某些高杂合或者特殊结构的样本上表现可能更好但代价是内存占用高得多且速度慢。只有在默认组装器拼出来的结果明显有问题比如线粒体被拆成很多短片段时才值得切过去试。内存是硬约束不是建议值。我见过最典型的失败场景是跑到组装阶段进程被系统直接杀掉日志里什么都没有只在系统日志里留下一行内存不足。这种情况你调多少线程都没用只能降数据量或者加内存。我常用的做法是把数据抽到组装器舒适区如果机器只有 32 GB 内存就把输入数据控制在 2 Gb 以内有 128 GB 就把上限提到 810 Gb。线粒体覆盖度早就够了多喂的数据只是在浪费时间。3.3 run_mode、dbnum、thread_number 的取舍这三个参数属于调了有感觉但不致命的类型理解它们能帮你优化时间不理解也不影响正确性。run_mode控制流程的彻底程度。快速模式会跳过一些耗时的搜索步骤适合先摸底、看这批数据有没有希望完整模式会做多轮钓取和回收结果更完整但慢。我的一般策略是先用快速模式跑一批样本看整体成功率对快速模式结果不理想的样本再单独用完整模式重跑。这样避免了所有样本都跑完整流程的时间浪费。dbnum控制把参考数据库切成几份并行搜索。这个参数跟线程数配合使用——数据库切得太少线程跑不满切得太多每个分片太小比对本身的固定开销占比上升。经验值是切成线程数的四分之一到二分之一比如 20 线程切成 510 份。thread_number就是线程数没什么好说的但有个细节这个值要同时满足组装器和比对程序的需求有些环境里它不能超过物理核心数超了反而变慢。3.4 两条可以直接抄走的命令说完参数给两条我自己常用的命令模板。第一条是完整的从头跑mitoz All \ --fq1 sample_R1.fq.gz \ --fq2 sample_R2.fq.gz \ --outprefix sample_mito \ --clade Chordata \ --geneticcode 2 \ --thread_number 20 \ --dbnum 8 \ --max_read_length 150 \ --run_mode 2 \ --fastp /opt/biosoft/fastp第二条是已经有序列、只想注释的情况——这是很实用的场景比如你自己手工拼好了环状序列或者从总基因组组装结果里抠出了一条 contigmitoz annotate \ --genefasta candidate_mito.fa \ --outprefix sample_mito \ --clade Arthropoda \ --geneticcode 5 \ --thread_number 10第二条命令省掉了组装环节几分钟就能出结果非常适合做参数试错。不确定 clade 或 geneticcode 填什么的时候先用第二条命令拿同一条序列反复试看哪组参数注释出的基因最全、最完整再去跑完整流程。这个技巧帮我节省了大量时间。4. 结果怎么读哪些文件有用哪些是噪音4.1 输出目录速查表一次完整运行会在输出前缀目录下生成一大堆文件第一次看会懵。实际需要关注的只有几个文件/目录内容重要程度前缀.fa最终组装的线粒体序列最关键前缀.gbGenBank 格式注释文件最关键前缀.summary.txt序列与注释的汇总统计高前缀.circos.pdf / .svg环形图中出图用前缀.gene.order.png基因顺序比较图中验证用前缀.depth.stat覆盖度统计高判断质量用前缀.trna.txt / .rrna.txt / .gene.txt各类基因的分项结果高核对用中间目录质控、组装、比对、注释、绘图各阶段中间文件低排查时才看我的习惯是先看.fa的长度再看.gb里基因数然后看覆盖度统计最后才看图。这个顺序能最快判断这条序列值不值得继续投入时间。中间目录不要急着删。排查问题时它们就是现场证据尤其是组装目录和比对目录能告诉你钓取这一步到底有没有命中以及命中了什么。4.2 判断组装质量的四个硬指标拿到.fa之后用四个指标快速打分基本能判断这条序列能不能用。指标一总长度是否落在该类群的正常区间。脊椎动物多在 16.5 kb 上下昆虫多在 1517 kb线虫偏短1314 kb软体动物变异最大从 14 kb 到 20 kb 都有。如果拼出来只有 8 kb那肯定是断了如果拼出 25 kb那可能混进了重复序列或者 numt。指标二覆盖度是否均匀。打开.depth.stat看沿序列的覆盖度分布。一条健康的线粒体序列覆盖度应该是平缓的高原偶尔在控制区有个小起伏控制区有串联重复比对时容易堆叠。如果出现一段覆盖度突然掉到个位数甚至 0说明那里有断点或者组装错误。指标三是否能成环。真正的线粒体基因组是环状的。判断方法是看序列两端是否存在重叠——如果末端几十到几百个碱基和开头高度相似那基本就是环状结构被线性化时的产物手动切掉重叠加起来就行。反过来如果两端完全对不上、覆盖度也不衔接那就是断的需要补洞。指标四基因数是否完整。动物线粒体基因组成非常保守13 个蛋白编码基因 22 个 tRNA 2 个 rRNA 37 个。少一两个可能是注释灵敏度问题比如某个 tRNA 反密码子特殊少五六个就是序列本身的问题了。4.3 注释结果的核对逻辑13 22 2数量对上了不代表质量对上了还得看细节。我核对的顺序是固定的先看蛋白编码基因。逐个检查三件事长度是否合理、起始密码子是否正常、有没有提前终止。有一个例外必须记住——有些基因在序列上确实会出现提前终止或者不完整末端因为它们的转录本要经过多聚腺苷化来补齐这是线粒体的正常现象不是错误。判断方法是把序列翻译出来看是不是只有一个位置异常、其余部分完全连贯。如果断得零零碎碎那就是真错了。再看 tRNA。22 个 tRNA 覆盖 20 种氨基酸其中亮氨酸和丝氨酸各有两个反密码子不同。少了某个 tRNA先别急着改序列可以先换一个 tRNA 识别工具交叉验证。线粒体 tRNA 二级结构短、变异大单一工具漏检很常见。低灵敏度模式或者换个算法往往就能捡回来。最后看 rRNA。两个 rRNA大亚基和小亚基的边界靠比对和二级结构确定自动注释经常有几十个碱基的偏差。这个偏差在系统发育分析里影响不大但如果要提交数据库最好手工修一下。还有一个必查项基因之间的重叠与间隔。动物线粒体里有几个固定的重叠模式比如相邻的两个 ATP 合酶亚基基因常重叠几个碱基、两个 NADH 脱氢酶亚基基因也常重叠。这些重叠如果在你的注释结果里消失了说明边界定错了。反过来如果出现大片不该有的基因间区也可能是边界问题或者拼装错误。5. 避坑实录报错、烂结果与补救办法5.1 运行期报错的排查顺序报错不可怕可怕的是没日志。我遇到过的报错按频率排下来大致是这几类排查顺序也有讲究。第一类找不到命令或者依赖。报错信息里会出现command not found之类的字样。这是环境问题八成是依赖没装全或者不在 PATH 里。解决办法是用参数显式指定路径或者把依赖目录加进环境变量。第二类内存不足被系统杀掉。特征是没有报错、进程莫名消失、退出码异常。查系统日志能看到内存耗尽的记录。对策是降数据量或者加内存别指望调参数能绕过去。第三类比对或搜索返回空结果。这类最麻烦因为流程可能正常结束但什么都没注释出来。原因可能是 clade 填错、参考库版本不对、序列本身不是线粒体。排查方法是去中间目录里翻比对输出文件看是不是真的一条命中都没有。第四类绘图阶段崩溃。前面全跑完最后一步挂掉。九成是绘图软件的环境问题Perl 模块缺失、R 包缺失。这类问题不影响数据本身.fa和.gb都是完整的可以手工用别的工具出图。心得养成先看中间目录再问人的习惯。一个跑失败的流程中间目录里留下的信息量远大于一句报错了。5.2 组装出来一堆 contig 怎么办这是最典型的结果不佳场景.fa里躺着七八条 contig长短不一加起来才勉强够线粒体长度。原因就那么几个我按可能性排序。原因一核基因组里的线粒体假基因numt干扰。这是头号杀手。动物的核基因组里普遍插入了线粒体片段的拷贝这些片段和真线粒体高度相似组装时会互相抢reads导致线粒体被拆成多个片段甚至拼出嵌合序列。怎么识别 numt三个特征覆盖度明显偏低numt 的覆盖度等于核基因组覆盖度通常比线粒体低一到两个数量级、含有提前终止密码子或移码因为插入后不再受纯化选择约束会积累突变、碱基组成偏离真线粒体通常 AT 含量很高numt 会向核基因组的 GC 水平漂移。用 reads 回贴到每条 contig 上看覆盖度一眼就能分出来。原因二数据量不够。特别是保存久的标本、非损伤取样如粪便、毛发的样本线粒体 reads 本来就少。这种情况下补测序是唯一解靠调参数没用。原因三控制区太难拼。控制区富含串联重复长度变异大是所有组装器的共同噩梦。这里我要说一句可能反直觉的话控制区拼不完整很多时候不影响使用。如果做的是系统发育或者基因组成分析把控制区标注为不完整区域即可不必死磕。原因四样本里混了别的物种。比如寄生虫、共生菌、食物残留。这类序列的 GC 含量和密码子偏好跟目标物种差得很远比对一下就能发现。对策是在组装前做去宿主或者按 GC/覆盖度筛 contig。补救的做法把目标 contig 挑出来用 reads 回贴看覆盖度把明显是 numt 的剔除然后对剩下的 contig 做延长拼接或者设计引物做 PCR 补洞。MitoZ 的注释模块可以单独对拼好的序列再跑一遍不用从头重来。5.3 注释缺基因、多基因、边界错位缺基因最常见的是缺 tRNA 和缺那个最短的 ATP 合酶亚基基因。前者是识别灵敏度问题换工具或者降低阈值能救后者是因为它实在太短很多类群只有五六十个氨基酸比对打分上不去被过滤掉了。对策是手工按位置和保守区找或者降低比对阈值重跑。多基因的情况通常是同一个基因被注释了两遍——典型原因是序列里存在重复片段numt 残留或者组装错误导致某段重复或者基因重叠区域被重复计算。看到 38 个、39 个基因先查重复。边界错位最隐蔽。表现为基因长度多了或少了十几到几十个碱基起始密码子不对。这种情况一般不影响下游分析但如果要提交数据库还是手工修一下比较稳妥。修的依据是翻译成蛋白后和近缘物种的同源蛋白比对看哪个位置开始出现连续的错配。还有一个特殊情况值得单独说基因重排。某些类群的线粒体基因顺序和标准顺序差别很大比如部分软体动物、部分寄生虫。如果注释结果里基因顺序乱得离谱先别怀疑工具去查一下这个类群是不是本来就有重排。这类信息在基因顺序比较图上能看得很清楚。5.4 常见问题速查表把上面这些浓缩成一张表出问题时按症状查症状最可能的原因处理办法流程跑完但注释出 0 个基因clade 填错或参考库不匹配换 clade 用注释模块单独重试蛋白编码基因普遍带提前终止geneticcode 填错按类群换密码子表重跑注释组装结果为多条短 contignumt 干扰或数据量不足按覆盖度筛 contig必要时补测序某段覆盖度骤降到接近 0组装断点用 reads 回贴定位设计引物补洞总长度远超正常范围混入重复序列或嵌合检查末端重叠切掉冗余部分tRNA 数量少于 22识别灵敏度不足换工具交叉验证降低阈值缺少最短的 ATP 合酶基因序列太短打分不够按保守区手工定位绘图阶段崩溃绘图环境缺失单独用其他工具出图不影响数据跑完进程莫名消失内存不足降数据量或加内存6. 收尾从自动结果到能提交的序列6.1 断点定位与补洞的实际做法自动流程给出的序列只能算半成品。要变成一条干净、完整、能提交的序列还需要几步手工活。第一步确认断点位置。把原始 reads 回贴到组装序列上看覆盖度在哪里断掉。断点通常表现为覆盖度急剧下降或者出现一堆 soft-clip 的 reads。把断点前后各一两百个碱基记下来这就是需要补的区域。第二步尝试延长。很多时候断点并不是真的没 data而是组装器在那里放弃了。可以把断点两端的序列当锚点去原始 reads 里捞跨越断点的 reads做局部延长。或者用另一套组装器重跑一遍两套结果互相补。第三步PCR 补洞。如果数据里确实没有跨越断点的 reads通常是重复区或者极端 AT 富集区那就只能回到实验台。在断点两侧设计一对引物扩增出产物做 Sanger 测序把序列接上去。这一步听起来麻烦但一条 16 kb 的线粒体通常也就一两个洞一两天能搞定。第四步确定起点和方向。数据库提交和后续分析通常要求统一的起点。常见约定是以某个特定基因或者 tRNA 作为序列开头。取哪个起点看你的研究类群惯例也可以参考近缘物种已提交的序列。定好之后把环状序列旋转一下切掉末端重叠得到一条标准的线性表示。提示补洞前的序列和补洞后的序列一定要都留着比对一下确认改动区域没有引入新的错误。我有一次补完洞忘了核对方向结果整条序列反向互补了白折腾半天。6.2 提交前的自查清单序列收拾干净之后提交之前我会过一遍这几项每一项都对应过真实的返工长度是否在类群正常区间是否确认为环状末端重叠已处理基因数是否为 13 22 2每个蛋白编码基因的起始、终止密码子是否合理是否存在异常的大段基因间区超过几百个碱基的空白需要解释序列的碱基组成是否符合该类群特征AT 含量、GC skew控制区是否标注清楚不完整是否说明注释文件的格式是否符合数据库要求是否有 numt 嫌疑片段残留其中异常大段基因间区这一项最容易被忽略。动物线粒体基因排布非常紧凑基因间区通常只有几十个碱基。如果你发现某两个基因之间空了两千个碱基那几乎可以肯定是要么拼错了要么混进了外源序列要么那段是没注释出来的基因。6.3 出图与二次美化自动生成的环形图够用但正式发表或者做汇报时往往需要再加工。我的经验是分工明确用自动工具出准确的图用绘图工具出好看的图。如果自动环形图能正常出直接拿矢量格式PDF 或 SVG进矢量编辑软件微调改字体、调颜色、补图例比重新画一遍快得多。如果自动出图失败可以退一步用在线注释工具重跑注释再出图或者自己用绘图代码画——把注释结果里每个基因的起止坐标读出来画成同心圆弧工作量不大且完全可控。基因顺序比较图的价值容易被低估。它不只是一张配图在判断组装和注释是否正确时非常有用把目标物种和几个近缘物种的基因顺序排在一起如果目标物种出现了大段倒位或者易位要么是真的有重排值得写进文章要么是拼装或注释出了问题需要回头查。两种情况都需要你重视。最后分享一个我自己常用的判断习惯拿到任何一条自动组装的线粒体序列先把 13 个蛋白编码基因翻译出来把氨基酸序列拼在一起看一遍。如果中间出现了莫名其妙的杂合、连续错配、或者某一段明显不像同源序列那问题就在那段。这个方法比看任何统计指标都直观几分钟就能定位到具体位置。另外工具本身在持续更新参数名和默认行为会变。我现在的做法是把每次跑通的命令和对应的版本号记在一个文本文件里下次直接复制。这个习惯帮我省掉过无数次上次明明能跑的困惑——因为上次跑的确实不是同一个版本。