只拿到一张组装指标表时,人很容易先挑 N50 最大、覆盖最全的一份。可是把两条染色体各剪开再互换一段,只要剪法对得巧,组装总长、序列条数、N50、参考覆盖率都可以一模一样。究竟哪项指标能识别这种变化?我们安装 QUAST,在同一份有两条染色体的参考上评估三套已知答案的组装。第一套全对,第二套交换两个跨染色体接头,第三套单纯打断一条序列。由此可以直接看到连得长、连得对、参考每一位有没有被覆盖这三种评价的区别。QUAST 5.3.0 手册
先把正确、错误接合和单纯断开分开设计
生成两条随机参考序列,甲十二万碱基,乙九万碱基。第一套组装是它们的准确拷贝,作为阳性对照。第二套将甲分成两个六万碱基的片段,乙分成三万与六万。把甲前半段接到乙后段形成一条十二万,再把乙前段接到甲后半段形成另一条九万。每一段都来自原始参考,碱基没有替换,也没有重复或遗漏,但两条序列各有一个跨染色体的人工接头。第三套只在甲的中间断开,分别保留两个六万及完整的九万乙,不增加任何连接。
将三套 FASTA 放进同一次 QUAST 5.3.0 运行,使用同一份已知参考、相同最短片段规则和比对参数。这样的设置可以避免一组因使用不同参考版本或不同长度过滤而显得更好。人工序列设计成无预设长重复,是为了让两处跨染色体连接能被清楚指认。在重复多、亲缘较远或存在真实重排的样本中,相同指标的生物学解释还得加上参考谱系和读段证据。这里的错误之所以可以确定,是因为我们自己知道每个碱基片段原来的连接位置。
为什么 N50 和覆盖率一起被“骗过”
N50 的计算只需要各序列长度。把片段从长到短排列,累计到组装长度一半时,当前片段的长度就是 N50。第一、二套的长度列表都是十二万和九万,总长二十一万,一半是十万五千。排第一的十二万已经越过门槛,所以两套都是十二万。中间到底接的是同一条染色体还是另一条染色体,长度列表里没有这项信息。第三套的列表是九万、六万、六万,第一条尚未过线,加到第二条时过线,N50 因而降为六万。
参考覆盖率问的是参考上有多少位置至少落入组装对齐中。三套都完整使用了原始参考的碱基,因此实际报告的参考覆盖率全部为百分之百,重复比也都为一点零零零。第二套乱接的两个片段分别仍能比回各自来源,它们的并集把参考覆盖得很完整。在一次批量组装里遇到百分之百覆盖同时还有重大错拼的提示,这两个数字是可以共存的。要验证错误到底在哪,需要查看每条查询的连续比对块及其相邻关系,不能期待覆盖率直接说出接头是否可靠。
| 组装 | 序列数 | N50 | QUAST 识别的重大错拼 | NGA50 | 参考覆盖 |
|---|---|---|---|---|---|
| 原样两条染色体 | 2 | 120000 | 0 | 120000 | 100% |
| 交换两处跨染色体接头 | 2 | 120000 | 2 | 60000 | 100% |
| 仅把甲分成两段 | 3 | 60000 | 0 | 60000 | 100% |

若比较多份来自不同亲缘背景的样本,先给参考序列建立同一版清单,注明是否包含替代序列、细胞器和未定位片段。只要改变参考总长或纳入的参考集合,达到一半所需的累计长度就可能改变。对短片段较多的组装,还要记录最短纳入阈值。这样读者才能按报告里的原始长度重新算出每个数,并判断两个样本的差异是否发生在同一评估范围。
实跑经验|把 NGA50 接回原始比对块
quast.py -t 2 --no-plots --no-html --no-icarus --no-gc -r reference.fa -o results intact.fa translocated.fa split.fa
awk -F '\t' '$1=="N50" || $1=="NGA50" || $1=="# misassemblies" || $1=="Genome fraction (%)"' results/report.tsv
QUAST 计算与参考比对有关的连续性时,会在相应的错拼边界拆分比对,再按对齐块长度统计。NGA50 用参考总长度作为达到一半时的分母。第一套保留了完整十二万,第二套两条人工连接被识别为两处重大错拼,拆开后 NGA50 只剩六万。第三套同样是六万,但并没有错误连接。于是即使拿到了一个更重视正确性的连续指标,仍要与错拼数、未比对长度一起看,才能区分“错误接长再被拆开”和“原本就很保守地留了断点”。QUAST 作者项目与参考评估工具
这次还有一个容易忽略的细节。第二套中最长的对齐块报告为六万零一个碱基,而人工拼接时制作的长块都是六万。局部比对可以在边界延伸,恰好相同的相邻碱基或比对边界选择会使报告区间与手工切点相差一位。没有逐位检查前,不能直接把这一个碱基解释成发生了一次插入。若要报告真实断点,应找出查询中的左右锚点、参考落点和方向,检查两侧是否有一小段可供多种等价放置的序列,并把可能范围保留下来。
实际使用时可先给每一处重大错拼提示分类,记录两侧落在同一染色体还是不同染色体、方向是否改变、间隔多大,以及原始长读段是否完整跨过。遇到同物种多个品系,软件报告的参考不一致可能正是研究对象携带的真实重排。一次原始跨接读段不足以消除所有歧义,最好检查独立分子的起点和重复背景,再决定该位置是组装改错,还是应当保留并作为结构差异报告。这样一套高 N50 的组装仍然值得用,但究竟哪里可以直接进入注释,哪里需要补证据,会有明确的坐标可以追踪。
下载本篇实验脚本与结果(12 篇合集 ZIP),解压后打开 075-quast 目录。
需要给多个样本或多版组装保存相同比对条件、断点表和原始 FASTA 时,可使用云生信的生物信息分析服务器批量复核。
