做完一轮组装抛光,原来两个怀疑有错的位点都得到了纠正。准备交付时,却发现总序列长度也变了。先别把差值都写成组装错误被修复。Racon 可以改写有读段支持的局部共识,默认行为还会影响哪些序列出现在输出文件里,以及窗口边界怎样修剪。这三个问题可以分别制造对照。我们造了两条长度已知的草图,只给其中一条的前部提供读段,并把三个人工错误埋在不同深度处。这样输出里的每个变化都能回到具体的原始碱基。Racon 项目说明
三份输入在生物学上必须对应起来
Racon 需要提供测序读段、读段对目标序列的比对,以及待抛光的目标序列。次序有实际含义。比对里的查询名称必须能在读段文件中找到,目标名称与坐标应当对应当前这版组装。多个项目都用类似的片段名时,比对文件容易用错,因此除了查名字,最好也查目标长度、参考文件校验值和上一轮修改日期。做第二轮抛光需要重新把读段比回第一轮的产物。只要某个窗口发生长度变化,原比对的后续坐标就可能平移,即使标题完全一样,也不能把上一轮的坐标直接挪用。
实验的第一条草图长六千四百个碱基,借助已知原序列,我们在第一千、二千五百和六千位各制造一个替换。第二条草图长一千六百个碱基,序列正确,但不给它任何读段。输入共有七十二条长度均为一千二百的无错读段,均来自第一条的前部。实际用 minimap2 比对后留下七十二条记录,最左起点为零,最右终点为四千七百零三,按左闭右开的坐标解释。前两处错误确实有读段跨过,最远的一处在全部读段的覆盖之外。这个安排将“碱基可能有错”与“现有观测能否纠正它”分开了。
把输出按名字、长度和碱基逐一清点
默认运行只输出第一条,长度从六千四百降为六千。我们按序列比对将输出重新放回原坐标,发现前两处替换都纠正了。未覆盖的第六千位依旧保留草图里的 A,已知真值为 C。长度差落在两处,一处是原序列开头的一百零三个碱基,另一处是紧接最后有读段覆盖的位置之后的二百九十七个碱基。两处相加四百,既解开了组装长度的变化,也解释了原来第六千位为何出现在输出的第五千六百位。把前后两行 FASTA 的同一列横着比较,会把一次裁切误读成它后面大量碱基全都错误。
第二轮在完全相同的读段、比对和草图上额外开启输出未抛光目标的选项。无任何读段覆盖的第二条回来了,长度仍为一千六百,与原序列逐字一致。第一条依旧是六千,那四百个碱基并没有一起回来。第三轮保留这个开关,再关闭窗口端部修剪。现在第一条恢复为六千四百,第二条也保留,前两处错误仍被改正,第三处仍为原来的错误字母。三个结果放在一起,可以分别追问是哪一个开关改变了序列成员,哪一个改变了第一条序列的长度。Racon 参数与窗口修剪说明
| 同一组数据的参数 | 有读段目标 | 完全无读段目标 | 原第 6000 位 |
|---|---|---|---|
| 默认 | 6000 bp | 无输出 | 仍为 A,输出第 5600 位 |
加 -u |
6000 bp | 1600 bp | 仍为 A,输出第 5600 位 |
加 -u --no-trimming |
6400 bp | 1600 bp | 仍为 A,输出第 6000 位 |

最省事的抽查位置有三个,分别挑在高覆盖内部、覆盖刚结束的边缘、整段没有读段处。在同一张记录里写出草图碱基、纠正后的碱基和独立读段支持,再另外登记参考坐标是否平移。这样即使共识文件没有逐碱基质量,也能指出每个判断依赖哪份证据。
实跑经验|如何给一次抛光做可解释的验收
minimap2 -t 1 -x map-ont -c draft.fa reads.fq > overlaps.paf
racon -t 1 reads.fq overlaps.paf draft.fa > default.fa
racon -t 1 -u reads.fq overlaps.paf draft.fa > keep.fa
racon -t 1 -u --no-trimming reads.fq overlaps.paf draft.fa > keep_full.fa
先导出前后两版的序列名和长度,分别查整条缺失、整条仍在但缩短、长度相同但碱基有改变三种情况。对第一种检查有没有任何读段对上,再查是否要求保留未抛光序列。对第二种先把两版相互比对,让差值落到原序列的某个区间,叠加读段覆盖的终点。对第三种才适合逐位检查替换是否得到读段或其他独立数据支持。记录深度时用实际对齐到该坐标的碱基计数,在窗口边缘还要查同一簇读段是否全都在附近同时结束。只有几条读段延伸到边缘,和内部多条读段完整跨过窗口,是两种不同的证据条件。
默认窗口长度为五百个碱基。本例后部差异恰好从最后一个仍被读段覆盖的区间边缘一直延伸到原序列第五千位,随后又保留了更远的原序列。这正是为什么不能把短少的四百个碱基统称为末端切掉四百。对含基因或结构注释的组装,要把原注释坐标迁移到新版序列后再核对,特别留意跨过剪切段的特征。使用保长选项也不能充当质量认证。没有数据覆盖的那处错误在三轮里都原样留下,保留它只说明输出继续携带原始草图信息。Racon 抛光流程实现
真实样本很少有从头到尾准确的真值。挑一部分独立获得的高准确读段或可信参考区间做复核,将有证据改善、证据矛盾和没有独立证据的区域分开记。遇到下一轮平均相似度略升,先确认比较的分母是否因为无读段片段退出而缩水。如果只在共同覆盖区间内比较,并另外公开退出了多少序列和碱基,读者就能判断改善来自纠错本身,还是来自评估范围变化。这个小实验也是设置停止迭代条件的一个实际起点。
需要替多份组装重复比对与抛光,并保留各轮序列和证据时,可使用云生信的生物信息分析服务器管理同一项目的多个版本。
