★ FEATURED ARTICLE 1. 动物线粒体基因组为什么值得单独跑一条流程线粒体基因组这块儿说大不大——动物 mtDNA 通常也就 15 到 17 kb一个环37 个基因看着像是随手就能拼出来的东西。但真把它塞进核基因组那套组装流程里跑结果往往很尴尬要么拼装器把线粒体 reads 当成低覆盖度噪声给扔了要么拼出一堆长短不一的 contig 压根不知道该信哪个要么干脆被核基因组里的线粒体假基因NUMT带偏最后交出去一条既不是线粒体、又不完全是核序列的缝合怪。我见过不止一次有人拿全基因组组装的默认参数跑动物样本最后在核组装里手动 grep 出一段 16 kb 的 scaffold 当线粒体用注释全靠外推审稿意见一来就顶不住。1.1 把 mtDNA 塞进核基因组流程里会出什么问题动物细胞里线粒体拷贝数极高肌肉、肝脏这类高耗能组织线粒体 DNA 占总 DNA 的比例能到百分之几血液样本可能连万分之一都不到。这个比例跨度决定了同一套流程换个样本就会翻车。核基因组组装器默认假设各区域覆盖度接近线粒体那条 16 kb 的环在这种假设下就是典型异常值覆盖度动辄高出核基因组几十上百倍k-mer 谱又和核序列高度相似因为核内存在大量同源片段。主流组装器要么把它拆碎要么把它嵌进核 scaffold 里当成重复区折叠掉。更麻烦的是核内线粒体插入片段。动物基因组里普遍存在把线粒体片段搬进核染色体的情况短则几百 bp长则几 kb偶尔还有近乎完整的拷贝。这些片段和真线粒体序列相似度能到 90% 以上但进化速率不同、常有移码和提前终止。一旦组装器抓到的是这些片段注释出来的蛋白编码基因就会出现内部终止密码子、长度对不上、tRNA 结构残缺。靠肉眼分辨极难必须靠覆盖度、比对分布和分类学过滤三件事交叉验证。所以线粒体基因组组装不是核基因组流程少给点数据它需要独立的一条链路先用轻量组装器从 reads 里把线粒体 contig 捞出来再用分类学信息过滤掉非目标序列最后用专门的注释引擎做基因边界识别。MitoZ 就是为这条链路设计的。1.2 MitoZ 在这条链路里到底替你做了哪些事MitoZ 把整条流程拆成组装和注释两块可以一条命令跑完也可以分开跑。组装端它默认用轻量组装器对 reads 做快速拼接拼接前会按你指定的类群做一次 reads 层面的过滤把明显不属于该分类群的序列剔掉拼接后再对 contig 做比对筛选留下能和目标类群线粒体参考序列对上的片段。这一步是它相对通用组装器最关键的区别——过滤发生在 reads 层面和 contig 层面两道而不是拼完再猜。注释端做三件事蛋白编码基因用比对加模型推断确定外显子边界tRNA 用专门的 tRNA 搜索程序加二级结构预测rRNA 用隐马尔可夫模型在结构域水平上找。三套结果最后合并成一个 GenBank 文件附带注释表、序列文件和环状图。它解决的痛点是标准化。动物线粒体的基因顺序、密码子表、tRNA 命名都有一套相对固定的规范手工注释十几个样本还能撑样本量上到几十上百没有一套可复现的流程基本没法交付。MitoZ 的价值不在于它的算法比人工精准而在于它把能交付的基线结果变成了流水线产物人工只需要在异常处做修正。1.3 什么时候不该用 MitoZ说实话MitoZ 也不是万能钥匙。三种情况我建议绕开或者换思路。第一种样本里线粒体 reads 占比低到千分之一以下。比如某些保存方式不当的组织、降解严重的样本、或者本身就混了大量外源 DNA 的环境样本。这种情况下组装端大概率只能吐出几 kb 的碎片后头注释再准也没用。得先想办法富集比如做线粒体分离、或者用长读长数据补。第二种研究对象是非动物类群。MitoZ 的参考库和密码子表默认面向后生动物真菌、植物的线粒体基因组结构差异太大植物 mtDNA 动辄几百 kb、还有大量重组硬套只会得到一堆错误注释。第三种需要精确到单碱基的群体遗传学分析。MitoZ 组装的 consensus 序列在多态位点上会做取舍如果研究对象是群体内的变异位点还是得回到 reads 层面做变异检测别拿组装序列当参考。提示MitoZ 输出的序列是 consensus不代表样本内所有线粒体拷贝的完全一致。动物个体内线粒体异质性虽然不高但在某些类群里确实存在做进化分析时要留个心眼。2. 装得起来和装得能用环境与参考库的两道坎MitoZ 的安装是新手最容易卡住的地方也是我见过最多人半途而废的环节。它的依赖不是普通的 Python 包那种pip install一把梭而是要在系统里备齐十几个外部可执行程序其中还有几个因为授权问题不能随包分发必须自己搞定。很多人装完跑第一条命令就报command not found然后开始怀疑人生。2.1 conda 装完之后仍旧缺的外部依赖最省事的装法是走 conda 频道主流生物信息学频道里有打好的包能一次性把组装器、比对工具、结构预测工具、序列处理工具都拉下来。装完之后先做一件事把可执行文件路径确认一遍。conda create -n mitoz_env python3.8 conda activate mitoz_env conda install -c bioconda mitoz which mitoz mitoz --help--help能出来说明主程序通了但这只是第一步。真正跑起来还会调用一串外部程序包括序列比对、组装、tRNA 二级结构预测、蛋白结构域搜索、多序列比对、格式转换等。有些装在 conda 环境里能自动找到有些需要手动放进 PATH。最典型的坑是基因模型推断程序。这个工具因为授权限制不能随 conda 包分发必须自己申请或从官方渠道拿到二进制文件解压后把可执行文件所在目录加进 PATH。export PATH/your/path/to/genewise/bin:$PATH genewise -h # 能打印帮助就说明通了第二个坑是环形图绘制。MitoZ 的环状图依赖 R 里的绘图包如果只装 Python 侧而不装 R 侧跑完注释后画图那一步会静默失败日志里只有一行不起眼的报错。要么提前把 R 和对应包装好要么干脆在命令里关掉画图回头用别的工具出图。第三个坑是多序列比对工具。注释端在确定蛋白编码基因边界时会调用几个比对程序做交叉验证缺一个可能导致部分基因注释不上而且不会明确报错只在结果里表现为基因数量偏少。我的建议是把所有依赖一次性验证一遍再上数据for cmd in mitoz megahit bowtie2 bwa blastn hmmsearch \ tRNAscan-SE arwen mafft muscle seqkit; do printf %-16s $cmd; command -v $cmd || echo 缺失 done这张表跑出来缺什么补什么比跑到一半再排查省事得多。2.2 本地参考库构建与磁盘规划MitoZ 注释端需要一个本地参考库用来做分类学过滤和同源比对。这个库不是装完就有的得自己下载分类学数据再构建。整个流程下来占用的磁盘空间不小解压后的分类学数据动辄十几 GB构建出来的库还要再占一份。构建的基本思路是先用官方脚本把分类学数据下载下来再用另一个脚本把它转成 MitoZ 能用的格式。python /path/to/MitoZ/annotate/download_taxonomy.py python /path/to/MitoZ/annotate/make_mito_db.py \ -i ./taxonomy \ -o ./MitoZ_db具体脚本路径和参数在不同版本之间略有差异建议直接去安装目录下的annotate子目录里看一眼文件名再执行。构建过程中会做大量序列提取和索引耗时从几十分钟到几个小时不等最好挂在后台跑。磁盘规划上我给的经验值是原始数据 中间文件 参考库每 100 个样本至少预留 200 GB。中间文件里最占地方的是比对索引和组装中间产物跑完一批后可以清掉。还有一点容易被忽略工作目录不要设在网络盘或自动同步的目录上。组装阶段会密集读写大量小文件网络盘的 IO 延迟会让整体时间翻好几倍同步盘还可能因为文件锁导致程序卡死。本地 SSD 是最稳妥的。2.3 clade 与 genetic_code 的选型逻辑这两个参数是整个流程里最容易设错、设错之后后果又最隐蔽的。程序不会报错只会给你一份看起来正常、其实基因全长对不上的注释结果。类群参数决定两件事reads 层面按哪一类群的参考序列做过滤contig 层面按哪一类群的数据库做比对筛选。选错了会出现两种典型表现过滤太严所有 reads 都被判为非目标组装结果为空过滤太松或选错大群混进来的近缘物种序列被当成目标留下。常见的取值包括后生动物整体、节肢动物、脊索动物、线虫、软体动物等。如果你的研究对象是昆虫老老实实选节肢动物比选后生动物整体过滤效果好得多因为后者的参考集太宽泛。分类层级越贴近研究对象越好但前提是这个层级在数据库里确实有足够的参考序列。遗传密码表决定蛋白编码基因的翻译规则。动物线粒体的密码子表和核基因组长得不一样而且不同类群之间还有分化。脊椎动物线粒体用一套多数无脊椎动物用另一套棘皮动物、扁形动物又各有各的规则。选错的直接后果是注释出来的蛋白编码基因内部出现一堆终止密码子基因长度明显偏短。类群常用密码子表编号典型异常表现脊椎动物2相对稳定选错时出现频繁内部终止多数无脊椎动物5选成 2 时丝氨酸/色氨酸编码错乱棘皮动物、部分扁形动物9选错时整条基因读不通部分寄生扁形动物14、21起始密码子识别失败实操上的判断方法很简单跑完之后看注释表里蛋白编码基因的长度如果多数基因比参考物种短一截或者内部终止密码子数量异常先怀疑密码子表。改一个参数重跑注释几分钟就能验证。注意不同版本的 MitoZ 对--clade的可选值不完全一致跑之前用mitoz all --help确认一遍当前版本支持哪些取值比照文档抄命令靠谱。3. 组装阶段从双端 reads 到候选 contig 的逐参数走查数据准备这一步看似没什么可说的但线粒体组装对数据量和数据质量的要求和核基因组完全是两个逻辑。核基因组是数据越多越好线粒体组装是数据量够用就行、但线粒体 reads 的绝对比例必须达标。这两者不能混为一谈。3.1 先算一笔账数据量到底够不够动物线粒体基因组按 16.6 kb 算目标覆盖度按 100 倍算需要 1.66 Mb 的线粒体序列。如果样本里线粒体 reads 占 clean data 的 0.2%那需要的总数据量是 1.66 Mb ÷ 0.2% ≈ 830 Mb。如果占比只有 0.05%就需要 3.3 Gb。这个账一定要在开跑之前算因为跑完发现数据不够返工成本很高。不同样本类型的线粒体占比差异极大新鲜肌肉组织能到 1% 到 5%肝脏类似全血通常低于 0.1%保存多年的标本、粪便、环境 DNA 样本可能只有万分之几。# 估算线粒体 reads 占比的快速方法 # 先取一部分 reads 比对到近缘物种的线粒体参考序列 seqtk sample -s 11 clean_1.fq.gz 200000 sub_1.fq seqtk sample -s 11 clean_2.fq.gz 200000 sub_2.fq bowtie2-build ref_mito.fa ref_mito bowtie2 -x ref_mito -1 sub_1.fq -2 sub_2.fq -p 8 \ | samtools view -bS - | samtools flagstat -比对上的 reads 数除以总 reads 数就是粗估的线粒体占比。这个方法在近缘物种参考序列可得的情况下很准如果连近缘参考都没有那就只能凭样本类型经验判断多给点数据总没错。数据质量预处理这一步我倾向于轻处理。接头必须去低质量末端必须切但不要做激进的长度过滤和纠错。线粒体基因组 AT 含量普遍偏高某些区域尤其是控制区本身就难测过度修剪会把这段序列整个切没最后组装出来缺一段。用主流预处理工具跑一遍默认参数然后把最短长度阈值设在 30 到 50 bp 之间别再往上调。3.2 一条完整的组装命令与参数含义先把数据跑一遍组装看看能拿到什么。mitoz assemble \ --fastq1 clean_1.fq.gz \ --fastq2 clean_2.fq.gz \ --out_prefix sample01 \ --workdir ./work01 \ --thread_number 16 \ --clade Arthropoda \ --genetic_code 5 \ --assembler megahit \ --insert_size 350 \ --requiring_taxa Arthropoda逐个说清楚--out_prefix决定所有输出文件的前缀也是结果目录的名字命名最好带上样本编号便于批量管理。--workdir是中间文件目录跑完可以整个清掉但建议先留着出问题时能从中间结果追溯。--thread_number按机器核数给一般给到物理核数的 80% 左右留点余量给系统。--clade前面说过选最贴近的类群。--genetic_code同上。--assembler指定组装引擎默认那个对线粒体这种小基因组已经够用换别的引擎可能耗时更久。--insert_size是文库插入片段长度从测序报告里查填错了会影响组装时 paired-end 关系的判断。--requiring_taxa是最终 contig 的分类学过滤条件和--clade不一定相同但通常保持一致。跑完之后先看结果目录里有什么。ls -lh sample01.result/ # 典型输出 # sample01.fa 组装序列 # sample01.gb GenBank 注释文件 # sample01.annotation.txt 注释明细表 # sample01.circos.pdf 环状图 # sample01.genes_info 基因位置信息3.3 拼出多段 contig 之后的收口手段现实里相当一部分样本第一次拼出来不是一条完整的环而是两到五条 contig。原因通常是线粒体基因组里有难拼接的区域AT 富集的控制区、串联重复、或者局部覆盖度极低的片段。先判断是不是同一条环被拆开。方法很简单把拼出来的 contig 互相做一次比对看看末端是否有重叠。makeblastdb -in sample01.fa -dbtype nucl -out contigs_db blastn -query sample01.fa -db contigs_db \ -outfmt 6 qseqid sseqid pident length qstart qend sstart send \ -evalue 1e-5 self_blast.txt如果 contig A 的末端和 contig B 的首部能对上且反向互补后也能对上说明它们大概率是相邻的。再看覆盖度bowtie2-build sample01.fa sample01_idx bowtie2 -x sample01_idx -1 clean_1.fq.gz -2 clean_2.fq.gz -p 16 \ | samtools sort -o sample01.bam - samtools depth -a sample01.bam | awk {s$3; if($30)n} END{print 平均深度,s/n}如果某条 contig 平均深度极高、另一条极低很可能是 NUMT 混进来了不是同一条环的两段。确认是同一环之后收口的做法有两条路。第一条是用 reads 在缺口区域做延伸把已有 contig 当参考把跨越缺口的 read pairs 提出来做局部组装。第二条更省事如果缺口很短几十到一两百 bp直接用 N 填充然后在注释阶段看是否影响基因边界。填充的写法# 按顺序拼接缺口处填 50 个 N cat c1.fa (printf gap\n; printf N%.0s {1..50}) c2.fa merged.fa更规范的做法是用专门的序列拼接工具把重叠区合并掉再手动检查连接处是否符合基因顺序。如果是蛋白编码基因被切开连接后必须能翻译出完整的开放阅读框这是最硬的验证标准。还有一种情况是拼出的序列是线性的、末端有重复。这说明环状基因组的起点落在了某个重复区组装器没法判断该在哪里收尾。处理方式是找重复序列把重复段手动截成一个剩下的部分首尾相接成环。3.4 reads 回贴与覆盖度均一性检查组装出环之后一定要做一次 reads 回贴验证。这一步很多人跳过结果交出去的序列在某个区域覆盖度为 0说明那段是组装错误或嵌合。bowtie2 -x mito_idx -1 clean_1.fq.gz -2 clean_2.fq.gz -p 16 \ | samtools sort -o mito.bam - samtools index mito.bam samtools depth -a mito.bam depth.txt # 分段统计找出低覆盖区 awk {binint($2/200); sum[bin]$3; cnt[bin]} END{for(b in sum) printf %d\t%.1f\n, b*200, sum[b]/cnt[b]} depth.txt \ | sort -n depth_bin.txt看这份分段覆盖度表理想状态下整条环的覆盖度应该相对均一波动在两三倍以内。控制区通常会有明显波动因为它 AT 含量极高、测序质量低、比对困难这是正常的。但如果非控制区出现覆盖度接近零的段落就要警惕。三类异常的对应判断覆盖度表现可能原因处理方式局部骤降为零组装错误、嵌合用 reads 局部重拼整条偏高数倍样本本身线粒体拷贝数高正常记录进质控表一条 contig 高、一条低低的那条可能是 NUMT单独做比对验证全程均匀但深度只有个位数数据量不足补测或放宽组装参数提示samtools depth默认会跳过比对质量低的 reads做覆盖度评估时按默认参数即可不需要额外放宽阈值否则会把错误比对也算进去。4. 注释阶段三类高频偏差的处理思路组装出一条完整的环只算走完一半。注释阶段的坑更隐蔽因为结果看起来像模像样但要真拿去投稿或者提交数据库问题就一个个冒出来了。我整理了三类最常见的偏差和对应的判断方法。4.1 蛋白编码基因的起始/终止密码子异常动物线粒体的蛋白编码基因有一套自己的规则起始密码子除了标准的 ATG还常见 ATT、ATC、ATA、GTG 等终止密码子除了 TAA、TAG还有不完整的 T 或 TA需要在转录后通过多聚腺苷酸化补全。注释程序如果严格按标准密码子表判断就会把这些基因判成不完整或者边界偏移。跑完注释先看明细表重点看两列基因长度和内部终止密码子数。awk -F\t NR1 $4CDS {print $1, $3-$21, $5} \ sample01.annotation.txt | head -20如果某个基因长度比同源物种的对应基因短了 20% 以上或者内部出现终止密码子要怀疑两件事密码子表选错或者基因边界被切错。密码子表的问题改参数重跑就能解决。边界的问题需要人工介入把目标基因的序列提出来和近缘物种的同源基因做比对看差异集中在哪个位置。seqkit grep -n -p cox1 sample01.fa cox1.fa blastx -query cox1.fa -db nr -outfmt 6 -max_target_seqs 5如果比对结果显示序列本身没问题只是边界差了几个碱基可以在 GenBank 文件里手动调整起始位置。MitoZ 输出的.gb文件是标准格式用文本编辑器就能改改完用格式校验工具检查一遍。4.2 tRNA 漏注释和补充策略动物线粒体基因组通常有 22 个 tRNA分布在两条链上其中丝氨酸 tRNA 有两个识别不同密码子亮氨酸 tRNA 也有两个。实际注释下来经常只能找到 20 个甚至更少漏掉的那几个通常是结构变异大、或者嵌在别的基因内部的。补注释的思路是双引擎交叉验证。MitoZ 默认调用了一个 tRNA 搜索程序但对结构异常的 tRNA 灵敏度有限。可以再用另一个基于二级结构预测的工具跑一遍把两次结果取并集再人工核对。# 用第二个工具独立扫描 tRNAscan-SE -A -o trna_scan.txt -f trna_ss.txt sample01.fa # 对比两个工具的结果位置落在这个并集之外的候选需要看二级结构图。动物线粒体 tRNA 有个特点很多缺少标准的 TψC 臂或 D 臂二级结构呈现三叶草退化形态。用通用 tRNA 模型去扫自然扫不到。这时候要放宽参数或者直接按基因间区的位置做推断——两个相邻蛋白编码基因之间的间隔区如果长度在 60 到 75 bp 之间大概率藏着一个 tRNA。漏注释本身不会导致分析做不下去但会影响两件事一是提交数据库时基因数量对不上审稿人会问二是做系统发育分析时少了一批可用位点。所以能补就补。4.3 rRNA、D-loop 与重复区的边界处理rRNA 基因在动物线粒体里就两个12S 和 16S位于两条链上。它们的边界相对固定但长度在不同类群之间有变化注释程序给的边界经常偏保守实际长度可能比输出结果长几十个碱基。判断方法是看 rRNA 两侧的 tRNA 位置。动物线粒体的基因排布高度保守12S rRNA 通常紧邻 tRNA-Phe16S rRNA 紧邻 tRNA-Val。如果注释出来的 rRNA 末端和相邻 tRNA 之间留了很大一段空白说明 rRNA 被截短了需要手动延伸边界到与 tRNA 相接。控制区也叫 D-loop 或 AT 富集区是整条环上最难处理的部分。它没有编码功能长度在类群间变化极大AT 含量能到 80% 以上还常含串联重复。注释程序一般只给一个大致的区域标记不给精确边界。处理原则是长度以能覆盖所有保守元件如中央保守序列、TAS 序列为准串联重复的数量按实际测序结果为准不同个体之间会有差异如果重复区导致组装长度不确定在 GenBank 文件里注明# 检查 AT 富集区位置 seqkit sliding -s 200 -W 500 sample01.fa \ | awk {gcgsub(/[GCgc]/,); print $1, 1-gc/length($0)} \ | sort -k2 -nr | head -5这段滑动窗口统计能快速定位 AT 含量最高的区域通常就是控制区所在。4.4 把结果整理成能交付的 GenBankMitoZ 直接输出的.gb文件能用但离能提交还差几步。需要检查的项目检查项合格标准常见问题序列长度与近缘物种相差在合理范围控制区长度异常基因数量13 个 PCG 22 个 tRNA 2 个 rRNAtRNA 缺失、rRNA 未标注基因顺序与近缘物种一致组装方向反向起始/终止密码子符合该类群规则边界偏移翻译结果无内部终止密码子密码子表选错基因顺序反向是个高频问题。动物线粒体基因组有固定的基因排列如果注释出来发现整条序列的基因顺序和参考物种正好相反说明组装出的序列是反向互补的。处理方法很简单seqkit seq -r -p sample01.fa sample01_rc.fa反转之后重新跑注释再核对基因顺序。选择起点也有讲究脊椎动物习惯以 tRNA-Phe 作为序列的第 1 位无脊椎动物各行其是参照同类群已发表基因组的惯例来定。注释完成后的可视化除了 MitoZ 自带的环状图我一般会再出一张线性基因箭头图方便在论文里展示。这类图用通用的基因组可视化工具都能画输入就是标准的 GenBank 文件。5. 四个踩坑场景的完整排查链路理论讲完说几个实际遇到过的情况。每个都按现象、假设、验证、结论的顺序还原方便你遇到类似问题时照着排查。5.1 只拼出 6 kb线粒体 reads 被稀释了现象一条命令跑完结果目录里只有一条 6 kb 的序列注释出来只有 5 个基因。第一假设数据量不够。先算账。原始数据 1.2 Gb clean data按公式反推需要至少 830 Mb假设占比 0.2%看起来够。但这是假设占比实测一下。验证抽样 20 万对 reads 比对到近缘参考比对率只有 0.03%。按这个比例算1.2 Gb 数据里只有 360 kb 的线粒体序列理论上只能覆盖 22 倍而且分布不均能拼出 6 kb 已经算运气好。第二假设样本本身线粒体含量低。查样本来源是保存了两年的乙醇固定组织降解严重且 DNA 提取效率低。处理两个方向。一是加大测序量按 0.03% 的占比推算要拿到 100 倍覆盖需要 5.5 Gb 数据成本不低。二是改变策略用多重 PCR 或捕获探针做线粒体富集把占比提到 1% 以上再测。实际选了后者富集后占比到 2.3%200 Mb 数据就拼出了完整环。教训数据量这个账必须在开跑前算而且要基于实测占比而不是经验值。低质量样本尤其要先做这个评估。5.2 拼出两条完整环NUMT 还是污染现象组装结果里出现两条长度都接近 16 kb 的序列分别注释出 30 多个基因都能闭合。看起来像是两条不同的线粒体基因组。第一假设样本混了个体。查实验记录是单个体提取排除。第二假设NUMT。核内插入的线粒体片段有时能拼出完整环尤其是插入片段比较新、还没积累太多突变的时候。验证三条线同时查。# 1. 覆盖度对比 samtools idxstats sample01.bam # 2. 与近缘物种参考的相似度 blastn -query contig_A.fa -db nt -outfmt 6 sseqid pident length -max_target_seqs 5 # 3. 检查是否有移码和提前终止 seqkit translate -f 5 contig_A.fa | grep \* | wc -l结果很清晰contig A 的平均覆盖度是 contig B 的 40 倍contig B 的蛋白编码基因里有 7 个内部终止密码子且 contig B 与参考序列的相似度略低、差异位点分布均匀符合 NUMT 的进化特征。处理保留 contig A丢弃 contig B。如果要做严谨的说明可以把 contig B 单独提交标注为核内线粒体插入片段。教训拼出两条环不一定是好事覆盖度差异和基因完整性是最快的判别指标。如果两条覆盖度接近、基因都完整那才是真的遇到杂合或污染需要回到 reads 层面分离。5.3 注释表里 PCG 一半带终止符现象组装没问题完整环但注释结果里 13 个蛋白编码基因有 7 个报告内部终止密码子基因长度普遍偏短。第一假设密码子表选错。样本是某种甲壳类默认参数选了脊椎动物用的那一套。改成无脊椎动物常用的那一套重跑报错减少但没完全消失。第二假设基因边界识别偏移。抓一个报错的基因单独看。# 提取基因序列并翻译 seqkit subseq -r 3400:4900 sample01.fa | seqkit translate -f 5翻译结果里前 30 个氨基酸完全正常之后突然出现终止再往后又接上了。这种模式说明不是密码子表的问题而是移码——起始位置偏了 1 到 2 个碱基导致整条读框错位。验证把同一个基因用近缘物种的同源蛋白做比对看哪个读框能对上。blastx -query gene.fa -db nr -outfmt 6 qstart qend pident -max_target_seqs 3结果显示从某个位置往后重新对齐得分更高说明确实存在偏移。处理手动调整起始位置往前或往后挪 1 到 3 个碱基重新翻译直到没有内部终止。MitoZ 输出的注释结果里记录了每个基因的坐标改完坐标后在 GenBank 文件里同步更新。教训内部终止密码子这个现象要分两种情况对待。全基因组普遍出现八成是密码子表问题个别基因出现八成是边界偏移。分开判断能省很多时间。5.4 环状图生成失败现象注释跑完了结果文件都在就是没生成环状图。第一假设R 侧依赖缺失。查日志果然有一行报错说找不到某个绘图库。处理装 R 和对应包。conda install -c conda-forge r-base r-circlize装完重跑画图那一步还是失败。第二假设中文字体缺失。项目里给物种名用了中文绘图时找不到对应字体。改成英文名重跑图出来了。教训出图这一步的报错经常被忽略因为它在整个流程的最后而结果文件已经生成很多人就不管了。建议把画图单独跑一次确认无误再进入下一步。6. 归档与交付一条流程该留下哪些东西跑完一批样本最怕的就是过几个月再回来完全不记得当时用了什么参数、删了哪些结果、为什么这个样本的注释被手动改过。所以流程的最后一环是归档规范。6.1 目录、命名与版本记录我给每个样本建一个独立目录里面按阶段分sample01/ ├── 00_raw/ 原始数据软链接即可 ├── 01_clean/ 预处理后数据 ├── 02_assemble/ 组装结果 ├── 03_annotate/ 注释结果 ├── 04_manual/ 人工修正记录 ├── 05_final/ 最终交付文件 └── run.log 完整命令与版本信息run.log里必须记录三件事软件版本、完整命令行、每个参数的选择理由。第三件事最容易被省略但恰恰是最有价值的——半年后你只看到--genetic_code 5不会记得为什么选 5 而不是 2。版本记录用一条命令就能搞定{ echo 环境信息 mitoz --version conda list | grep -E mitoz|megahit|bowtie2|blast echo 运行命令 history | grep mitoz | tail -5 } run.log6.2 一张能自证的质控指标表交付时附一张表比口头解释一百句都管用。我通常记录这些指标指标记录值参考范围clean data 总量3.2 Gb视样本而定线粒体 reads 占比0.85%越高越好组装序列总长16,542 bp15–17 kb脊椎动物contig 数1环状1–3平均覆盖度168× 50×覆盖度均一性变异系数 0.31 0.5蛋白编码基因数1313tRNA 数2222rRNA 数22内部终止密码子数00有了这张表任何一个内行扫一眼就知道这条序列能不能用。审稿人问你怎么保证组装质量直接甩表。6.3 审稿人最爱问的那几个问题投出去之后被问得最多的三个问题提前准备好答案能省一轮返修。第一个是控制区长度为什么和参考物种差这么多。答案通常是串联重复区拷贝数不同这是真实存在的种内变异不是组装错误。应对方式是给出重复区的序列结构说明拷贝数范围。第二个是基因顺序为什么和已发表的不同。除了方向选反这种低级错误还有可能是类群内部真实的基因重排。动物线粒体确实存在基因顺序变化尤其是 tRNA 的位置。应对方式是把近缘几个物种的基因顺序列出来对比说明这是已知的重排类型而不是注释错误。第三个是NUMT 怎么排除的。这个必须给出数据支撑覆盖度对比、翻译完整性、与参考序列的相似度分布三样至少给两样。我把这些整理成一份随交付文件走的说明文档每次提交前对照检查一遍比临时抱佛脚强。说个我自己的体会线粒体基因组这条流程技术上真正难的地方不多难在细节的一致性和可复现性。同样一批样本A 样本拼出完整环、B 样本只有三块碎片、C 样本冒出两条环这种不一致才是日常。所以别指望一套参数跑天下真正有价值的是那套排查思路——先看覆盖度、再看分类学归属、最后看翻译结果三层筛下来绝大多数异常都能定位到具体环节。另外手动改过的每一条注释都要留痕迹哪怕只是在注释表里加一行备注半年后你会感谢当时的自己。
阅读完成 · 觉得有帮助？
