序列对比 python 序列对比算法

• 
  

        今天给大家介绍的是2016年在Nature Communations 发表的一篇三代序列比对算法GraphMap，基于gap-gram查找和k-mer图来进行查找和比对，思想很巧妙，对错误率有很好的的鲁棒性。因此，本次简单介绍GraphMap算法的主要思想，希望对大家有所启发，可以用GraphMap算法来处理自己的序列。

GraphMap主要步骤

        GraphMap主要分为四个步骤，如下图所示。整体来看，GraphMap同样遵循“种子-扩张”比对思想。即首先构建基因组的种子索引，然后基于gap-gram提取待比对序列的种子，构建k-mer图，找到一组线性种子链（chain），再对种子链进行首尾优化，最后对非种子区域进行详细序列比对，与种子区域合并组成完整的序列比对结果。

        其中索引阶段采用成熟的哈希或者BWT-FM进行索引，这个是目前绝大多数比对算法采用的技术，都可以快速的进行子片段（种子）查找。其最大的特点是采用采用gap-gram进行索引，即含有空位的k-mer，主要是考虑测序误差的存在。这一步其他方法都采用 k-mer 进行精确查找。而gap-gram允许空位存在，可以找到模糊索引，降低测序错误带来的影响。找到相似 k-mer之后，构建 k-mer图，从图的角度来寻找chains。与他方法相比，每一步都有自己的思想，但目的都是一样的，就是采取的策略有所不同。

构建基因组索引

        文章并没有详细介绍索引这部分，因为不是文章的主要内容，只简单提了句：当 k<10 的时候，用哈希索引，否则用后缀数组（Suffix Array）。在我看来，其实只要是速度不慢的索引技术，都可以拿来用。

gap-gram选取与查找

        根据相关文献介绍，Gapped q-gram是一些长度为 s (s > q)的模式，模式中除 q 个固定位置的字符外，其他位置的字符被认为是空位，称为(s, q)-gram。其中 s是截取的总长度，然后提取 s 中的 q 个字符，在匹配过程中只需匹配这 q 个固定的字符，而不关心空位上的字符。在我看来，其实就是变相的删除固定位置上字符的 k-mer。如下图所示。

        由上图可知，其中构建基因组索引的时候，先提取位置1-13的k-mer，然后把第七个碱基删除，余下12个碱基组成新的种子（seed），构建查询列表。

        对于待比对序列，如下图所示，同样提取三种模式的 k-mer，第一种和基因组一样，把第七个碱基删除掉。第二种将第七个和第八个碱基删除掉，第三种是把第13和第14个碱基删除掉。将提取得到的三种种子，通过之前构建的基因组索引表就可以查询在基因组中出现的位置。这样通过这三种模式有选择性的删除一些碱基，可以进行模糊查找而不是精确查找，避免了因测序错误而查不到匹配 k-mer的情况。

候选窗口选取

        这个方法与rHAT方法（rHAT，国内首个三代序列比对算法）类似，将基因组分割成一系列窗口，如下图所示，记录每个窗口查询到的gap-gram个数，然后选取个数比较多的窗口作为候选比对窗口。

构建 k-mer 映射图

        截取候选区域，构建k-mer映射图，如下图所示，每一个节点只与前三个节点相连。即每一个节点有三条（l -1条）边。其中 l 是考虑了测序误差的存在。构建好参考序列的 k-mer映射图后，提取待比对序列的 k-mer，把出现的 k-mer变成另一种颜色，然后连接在一起的 k-mer就可以组成待比对的框架，将成对待比对的序列分割成小碎片（又称非种子区域），只需比对小碎片，再与种子合并，最后组成完成的序列比对结果。 

 优化比对框架

        上一步得到的比对框架构建好后，再对其进行最后一步优化，如下图所示。其中左上图为候选区域和待比对序列的种子位置信息图，右上角是通过上一步得到的比对框架，然后文章用一个LCSk（Longest Common Subsequence in k Length Substrings）进行优化，选取更好的比对框架（左下角图），当然优化后的框架和没优化前的框架大体差不多。然后对非种子区域进行序列比对，与种子区域合并，最终得到详细的序列比对结果。

比对结果

        下图是GraphMap算法在三组模拟数据集上的比对结果，与其他五种方法进行比较：BLASR（三代测序序列比对利器-BLASR，更小更快更方便）、DALIGNER、BWA-MEM和LAST。如下图所示，可以看出GraphMap方法在比对碱精度和召回率均表现出良好的效果。希望可以对序列分析者有所帮助。剩下详细的比对结果大家可以参考其原文。 

代码下载地址：

https://github.com/isovic/graphmap

原文链接：

https://www.nature.com/articles/ncomms11307

参考资料：

1. Sović, I., Šikić, M., Wilm, A. et al. Fast and sensitive mapping of nanopore sequencing reads with GraphMap. Nat Commun 7, 11307 (2016). https://doi.org/10.1038/ncomms11307