Merqury 实跑｜同一份组装，为什么换一批评估读段，QV 从 34 降到 23？

假设你拿到同一个样本的两份质控报告，组装文件的校验和一样，一份报 QV 34，另一份只有 23。回头检查抛光程序，很可能查错了地方。Merqury 用准确读段为组装中的短序列作证。它既会察觉写错的碱基，也会察觉读段库没能替一段正确序列作证。我们在一条已知真值的序列上分别制造这两种情况，追踪总分、短序列计数和错误区间如何一起变化。Merqury 方法论文

先问哪个集合在做分母

把一条序列上连续的二十一个碱基当作一个词，窗口每次前移一位。和它对应的反向互补形式合并处理，这样从哪条链测到都能提供证据。组装里的词在准确读段库里也出现过，说明测序数据支持这个局部顺序。如果只在组装里找到，脚本把它计入缺少支持的数量。读表时先写清比较的方向。反过来查读段里有、组装里找不到的词，是追踪缺失序列的另一道题，分母、可信词筛选和生物学含义都变了。作者的 QV 计算脚本

一个孤立的碱基写错，会破坏盖住它的一串窗口。只要离两端足够远，长度为二十一的窗口最多有二十一个盖住这个位置。把二十一个异常词除以组装长度，会把同一处改变算上很多遍。脚本先估计一个完整窗口被支持的比例，再按窗口长度折回到单碱基错误概率，最后转成对数分数。相邻错误共享窗口，重复序列可能在别处获得支持，真实基因组很少满足最简单的独立模型。把分数说明为基于读段支持的估计，读者才能判断它能回答多细的问题。

实跑经验｜组装里始终只有一个替换

我们用固定随机种子生成长两千六百碱基的真值，复制一份，在从一开始计数的第一千三百零一位改写一个碱基。评估读段始终取自原始真值，每条一百五十碱基。第一组以五十碱基步长平铺到末端，共五十条。第二组分别铺在预先选定的缺口两侧，共四十四条，中间两百个碱基没有任何读段经过。两次运行读入同一份组装，均用 meryl 1.4.2 建库，再调用 Merqury 1.4.1 自带脚本。这个输入刻意消除了测序噪声，让缺少取样这件事可以单独观察。

这个特制序列的组装窗口恰好全部不同，因此分母两千五百八十也可以从序列长度减去二十直接复核。第一组只有跨过替换的二十一个词匹配失败，模型给出的碱基错误率约为万分之三点九，和真值中的一个错误除以两千六百相近。第二组多出二百二十个组装独有词，原来的替换仍在那里，但新异常来自抽样缺口。从上游再做一轮抛光，也无法凭空补回读段库里缺少的证据。

区间的边界尤其容易误读。第一组的异常区间长四十一个碱基，知道答案的我们却只改过正中一个位置。脚本报告的是连续异常词覆盖范围的并集，最左边的词从替换点上游二十位开始，最右边的词又向下游延伸二十位。第二组新增区间长二百四十个碱基，而人工留下的缺口只有二百个。把这份区间清单直接当作需要剪除的序列，会在两边一并丢掉有证据的碱基。应先回到原始比对或深度轨迹，查看内部与边缘分别发生了什么。

怎样复算报告上的数

meryl count threads=2 memory=1 k=21 output reads.meryl reads.fa
bash "$MERQURY/eval/qv.sh" "$PWD/reads.meryl" "$PWD/draft.fa" report
cat report.qv
cat report.draft.qv
cat report_only.bed

我们另外写了一个短程序，对正向词与反向互补词取字典序较小者，建立两个集合，直接数出只在组装里出现的元素。两种输入的数量都和脚本相等。再用组装独有词数除以总词数，取其补数的二十一次方根，用一减去这个根就得到估计错误率。以十为底取对数后乘负十，两组数字都能复算到小数点后四位。真实项目包含重复词时，应照所用版本脚本核对计数定义，不能把这条随机序列上碰巧相等的“窗口数”和“不同词数”直接互换。meryl 计数工具

怎样区分证据薄弱和组装错误

如果新一版组装总分下降，先固定样本、评估读段、过滤规则和词长，用同一个数据库重算新旧两版。把异常区间按染色体和局部深度整理，先看新出现且原始读段覆盖充足的部分。这些位置值得逐个查看读段支持的碱基、附近插入缺失和重复拷贝。若异常集中在多个样本都难测的片段，先查建库与测序的覆盖偏差。有独立文库时，在同一位置比较支持情况，比只把两份总分排个名次更能缩小排查范围。

二倍体组装还要核对评估数据来自哪个个体。用不同个体的高质量读段评价同一条单倍型，真实的个体差异也会留下组装独有词。反过来，一个拷贝被错误折叠，留下的那份序列仍可能处处得到读段支持。需要找这类问题时，可以同时检查可信读段词的完整性和拷贝数谱，并对可疑区间核实深度。对于缺少独立证据的地方，在记录中留下待验证区间及依据，通常比把整条染色体判成低质量更有用。Merqury 官方使用说明

实际交付时，我会把每条序列的表和全组装分数一起保存，并附上建库输入的校验值、词长和窗口计数。这样下游若问为什么某条短质粒拉低了分数，可以直接回到那一行，计算它的独有词对总数贡献了多少。两次评估若采用了不同的读段子集，报告应并排给出每套数据库的覆盖情况，再决定分数差异能支持哪一种判断。

下载本篇实验脚本与结果（12 篇合集 ZIP），解压后打开 074-merqury 目录。

需要在同一套准确读段上反复比较多版大型组装时，可以在云生信的生物信息分析服务器保留建好的词库，供后续版本直接复用。