MECAT: fast mapping, error correction, and de novo assembly for single-molecule sequencing reads

一、 简介

MECAT由中山大学的Chuan-Le Xiao团队研发，其论文于2017年发表于《nature》。

MECAT是一个基于比对纠错的组装工具,获取基因序列后将其分割为多个长读段，并建立索引表，通过其原创的比对算法获取读段间的overlapping，使用比对结果进行纠错操作，最终将纠错完成的长读段进行组装，输出一整条基因序列。

论文中提到，工具的主要创新点是： 提出一套伪线性对齐评分算法（pseudolinear alignment scoring algorithm）来过滤过多的对齐，节省计算资源,大大减少运行速度。

Chuan-Le Xiao团队针对MECAT的各项性能指标设计了大量的实验，最终结果清晰的说明：MECAT工具的比对速度数倍高于其他方法，敏感度与运行效率平衡极佳，所需运行资源非常少，非常适合进行规模较小的基因序列比对、纠错、组装工作。

论文地址

项目地址

论文补充材料地址

二、方法介绍

(一) Indexing and matching of reads 建立索引&寻找匹配

MECAT寻找读段之间的潜在匹配是基于读段之中的k-mers（长度为k的子字符串），长度为L的读段总共有L-k+1个k-mer。算法首先使用k-mers作为key，k-mer在读段片段中的位置作为value，建立哈希表，对整个读段进行索引。关于片段，算法将每个读段分割为多个block，每个block的长度为B，其值通常为1000~2000bp。

当我们有一条待比对长读段，需要找到与它存在匹配的(相同)读段所在的位置，就扫描待比对长读段中的k-mers，并用其于哈希表中获取匹配的k-mer的位置，这些k-mer位于其他区域。为了减少运行时间，算法对于待比对长读段中的k-mers进行采样：使用长度为sl的滑动窗口，沿着当前所在的block滑动，每次滑动时取从滑动窗口头部开始的k-mer，并将滑动窗口从头部移动到尾部。经过采样，最终得到的k-mers数量约为整个读段中k-mers总数的1/sl（sl的默认值为10）。如若一个待比对的block和另一个block之间存在匹配k-mers，且数量大于预定义的阈值m，则两个block被判定为相互匹配；若两个读段中至少存在一对block相互匹配，则这一对读段被判定为相互匹配。

论文中还论述了一对匹配的block中相互匹配的k-mers的期望数量，具体内容请自行阅读论文。

图1

(二) Filtering false-matched reads using the distance difference factor score 过滤过多匹配

MECAT提出了一种原创的伪线性评分算法来过滤大量的弱信息匹配读段，被过滤的block不进行计算考察。评分算法包含两个步骤：相互评分（mutual scoring）和拓展计算（extension scoring）。

在相互评分中，算法首先在相互匹配的读段中随机选择一对匹配的block，并block进行标记，然后对其中匹配的k-mers对进行打分。打分的公式如下：

\[DDF_{i,j}=|1-\frac{p_i-p_j}{p_i^{'}+p_j^{'}}|\]

其中$p_i,p_j$代表block中第i、j个k-mer的位置，而${p_i}^{‘}，{p_j}^{‘}$代表另一个block中第i、j个k-mer的位置。公式中分数代表了这一对k-mer在两个不同的block中相对位置是否相同/近似，若近似则可以认为这一对k-mer更有可信度。

对于一对k-mers的DDF值若小于阈值ε（默认设为0.3），则表示这两个k-mers相互支持对方的可信度，于是将两个k-mers的分数都加一，在所有k-mers都两两成对完成DDF计算后，将分数最高的k-mer作为这一block内的seed，以进行下一步操作。需要注意的是，若一个k-mers在哈希表中能找到多个匹配项，则不置入计算。

图 4

当一个block中出现了seed后，算法将开始对其进行拓展计算：将与该blcok相邻的的存在匹配的block1中的每一个存在匹配的k-mer，都与block的seed两两成对，计算DDF值，如若值小于阈值ε，则将seed的分数加一。当block1中至少80%存在匹配的k-mer与seed成对的DDF都小于阈值ε，则标记block1，此后不再对其进行评分计算。当一次相互评分与拓展计算循环结束后，读段中仍然存在尚未标记的存在匹配的block，则在其中随机选择一个开始评分计算。

图 5

相互评分与拓展计算的时间复杂度分别为$O(N^2)$、$O(N)$，其中N代表与k-mer的匹配数。由于相互评分过程中k-mers的数量较少，所以整个评分过程可以在伪线性时间复杂度内完成。

(三) Pairwise alignment of single-molecule sequencing reads. 单因子序列读段的成对比对

在成对比对中，算法将每个block的长度设置为2000bp。当两个读段完成了计算，所有block都进行了标记后，算法将读段中的k-mers按照分数进行排序，使用top-rank的mers作为seed，对这一对读段进行局部比对（local alignment）。如若两个单分子序列读段的重叠区域长度都大于2000bp，且重叠序列的错配率（mismatch rate）两倍小于读段的错误率（error rate），则判定这一对读段相互匹配，最后将比对结果输出。

(四)Aligning single-molecule sequencing reads to a reference genome. 单因子序列读段比对到参考基因组

读段比对到参考基因组的工作与读段两两比对较为类似，区别主要在于前者多了一个根据不同情况设置block长度的策略，以及一个针对k-mers进行采样的滑动窗口策略。

算法首先将参考基因组分成长度为B的block，并使用其中k-mers建成索引表，然后将读段也分成长度为B的block，对其中的k-mers进行采样，并在索引表中搜索这些k-mers，找出匹配的block，对其施加前面提到的比对策略。完成读段中所有block的评分操作后，将所有k-mers按评分进行排序，使用top-rank mers作为种子进行局部比对。

读段与参考基因组的比对工作有两个步骤以应对参考基因组过长导致的比对效果不佳等问题。在第一步中，算法将block长度B设为1000bp，k-mers采样滑动窗口的长度sl设置为20而进行比对，但有一部分block中存在匹配的k-mers数量较少，或是分布不均匀，所以无法找到匹配项。在第二步中，算法将B设为2000bp，将sl减小到10，然后重新比对在第一部中没有完成匹配的block。

由于在上述的两个步骤中，第一个步骤能找出读段内大部分block的匹配项，且第一步所需的计算资源小于第二步，所以这个策略能在保持比对精度的基础上减少资源占用。

该部分论文还对读段中局部结构的修剪策略做了论述，如有需要请自行查阅。

(五)Correcting single-molecule sequencing reads.单分子序列读段纠错

在MECAT中，读段纠正包含两个步骤，第一个是在读段中采用成对比对，第二个是根据一个读段的共识表来构建正确读段。下面进行具体说明。

在第一步中，算法把将要进行校正的读段称为template read。对于要进行校正的template read，将其与所有匹配读段进行DDF计算，把这些匹配读段根据其中的k-mers分数进行排序，然后从高到低进行局部比对。需要注意的是，在寻找匹配项的过程中，如果两个读段重叠部分的长度小于其中较短读段的90%，则舍弃这个匹配读段，以消除嵌合读段和重复子序列（chimeric reads and repeat subsequences）的影响。

文中提到，一旦收集到100个overlap或template read的所有匹配读段，比对操作就会停止，且由于DDF得分能粗略估计overlap的长度，所以根据DDF分数从高到底进行局部比对能提高效率，避免计算无信息的重复overlap，加快校正。关于DDF分数能够估计overlap长度的说法，我个人的理解是，这里的DDF分数指的是读段内所有k-mers的在比对操作过程中不断累加，最终得到的分数，分数由k-mers相互比对得来，值越高就说明k-mers数量越多，那么重叠长度也就越长，所以DDF分数才能粗略估计重叠长度。

图 8

在第二步中，MECAT结合了其他两个工具——DAGCon和FalconSense的原理，提出了一种自适应单分子读取错误校正策略。算法使用在第一步中template read与匹配读段进行局部比对所得到的信息构建了一个共识表（consensus table），其中包括读段每个位置上碱基的比对信息，这些信息是匹配、插入和删除操作的计数统计值。有上图中比对信息构成的共识表如下图。其中mat代表match，del代表deletion，ins代表insertion，base代表碱基。

图 6

在完成共识表后，MECAT根据其中的信息将template read分成三个区域：

不重要区域：$\frac{match}{match+deletion} > 0.8$ && $insertion < 6 $

删除区域：$\frac{deletion}{match+deletion} > 0.8$ && $insertion < 6 $

复杂区域： $insertion \geq 6$

通过分区获得的共识结果（consensus result）如下图所示。

图 8

算法针对三种区域有不同的处理方法，其中匹配区域完整保留，不做处理；删除区域直接去除，不留gap；复杂区域则使用DAGCon的原理，通过建立DOG来得到共识序列（consensus sequence）。获得共识序列的操作如下图所示。

图 7

根据共识区域将复杂区域填充后，得到最终序列，如下图所示

图 9

该部分论文还说明了MECAT优化硬件的策略，有需要请自行查阅。

(六)De novo assembly using single-molecule sequencing reads.从头组装单分子序列读段

MECAT组装策略分为三步：

  • 对于每一条read，将其与其他read进行成对比对，在这过程中不进行局部比对，而是根据DDF分数高低来选择出前100条与其最匹配的read
  • 对于每一条read，使用其匹配的reads来对它进行纠错
  • 对经过的纠错的reads进行成对比对，并将比对结果输入到 Canu (v1.0) 的Unitig Construction模块来将读段组装为unitigs.

MECAT工具安装较为繁杂，且其中存在BUG，建议直接安装MECAT2，集成度比较高，更容易使用。

MECAT2项目地址：https://github.com/xiaochuanle/MECAT2

谢谢，再会！