minimap2论文算法详解(主要针对RNA-seq)

Heng Li, Minimap2: pairwise alignment for nucleotide sequences,Bioinformatics, Volume 34, Issue 18, 15 September 2018, Pages 3094–3100, https://doi.org/10.1093/bioinformatics/bty191

前言

最近研究将三代测序RNA-seq比对到参考基因组，而Minimap2是三代测序数据重要的比对工具之一，因此想要了解minimap2所使用的比对算法，对其论文进行了研究。在读论文的过程中，发现很多公式参数不容易理解，也没有具体实例加深理解，网络上对算法的讲解也很少，在此记录一下自己对其中一些算法的理解，便于学习交流。

---

本文按照自己理解对论文进行解读，欢迎批评指正！

1. 引言

对研究问题背景及意义的一些简单介绍

1. 单分子实时测序技术（SMRT）和牛津纳米孔技术（ONT）产生读段长度超过10kbp，错误率为15%。最近测序技术的进入有望实现平均100kb的超长读段，高通量的全长mRNA或cDNA读段和基因组的contigs长度超过100Mb。——现有校准工具不能处理此类大规模数据或处理效率低下，这推动了新的校准算法的发展。
2. SMRT和ONT均已应用于剪接mRNA（RNA-seq）的测序,传统mRNA对准器（Spaln和GMAP）没有为长噪声序列读段进行优化，且比专门的长读段对准器慢几十倍。——Minimap2是第一款专门为长噪声读段设计的RNA-seq对准器，还扩展了原始算法，使其能比几种主流的短读段映射器更快。

——既对基因组DNA进行比对，也可以对mRNA进行比对。

2. 方法

1. 根据参考基因组建立minimizer索引表：将基因组序列的minimizer存储在哈希表中（minimizer指一段序列内最小哈希值的种子，因为所研究问题是将RNA-seq比对到参考基因组，因此，是对参考基因组建立索引）；
2. 对于每一条待比对序列（也就是三代测序的reads读段），找到待比对序列所有的minimizer，根据minimizer[x-w+1,x]中x的位置建立anchor表；
3. 根据最优chain分数，将anchor连接起来，从而在参考基因组上进行定位；
4. 再对非种子区域利用动态规划进行扩增。
   [注：这里参考基因组建立了一个minimizer索引表，每条read通过与minimizer的比对都有一个排好序anchor表，一个minimizer可能在reads找到很多匹配段，所以再根据打分，将符合要求的anchor连接起来形成chain，从而定位到参考基因组上。以上很多名词后面都会进行一一解释和计算]
   [图片仅供参考，感觉不是很准确]

主要流程如下：

2.1 计算minimizer

公式太难打了，我把PPT截个图吧hhhh

符号说明：

针对每个k-mer的哈希值的计算：

伪代码 by[Li,H. (2016) Minimap and miniasm: fast mapping and de novo assembly for
noisy long sequences. Bioinformatics, 32, 2103–2110.]：
举一个具体的例子来计算minimizer：

更多具体细节参考[Roberts,M. et al. (2004) Reducing storage requirements for biological
sequence comparison. Bioinformatics, 20, 3363–3369.]，放几张论文里的图加深理解：

最终，哈希表里存储的位置，是i+j，i就是w窗口所在位置，j就是在每个w块里minimizer的位置，所以最终的位置就是这个k-mer的minimizer在整个参考基因组中的位置。

2.2 chaining

主要流程：

2.2.1 找到最优chaining scores

算法如下：

将anchor定位到参考序列中，整个动态规划过程就是寻找符合连接标准的anchor，如果两个anchor位置比较近，也就是两个anchor之间有匹配的碱基，则对其进行加分；如果两个anchor之间存在gap，则对其进行罚分，通过这种方式可以找到相近的anchor，并将其连接起来，从而可以找到read每个anchor段在参考基因组上的位置。

2.2.2 回溯

算法如下：

通过动态规划算法对每个anchor进行打分，找到anchor之间的关系，然后进行回溯，将anchor连接起来成为chain。

2.2.3 识别primary chains

算法如下：

chain之间可能会有很多重叠，因此通过上述算法找到一个主串，且通过评估映射质量公式来评估primary chain的质量，chain质量越高越有利于得到好的比对结果。

2.2.4 评估每个序列的散度

算法如下：

就本人理解，此处序列散度是为了评估chain质量，散度越大越好还是越小越好。理解不是很透彻。

2.2.5 均聚物压缩k-mers进行indexing

算法如下：

为了节省存储空间的一种方法。

2.3 对齐基因组DNA

流程如下（因为RNA-seq在第一个上面改的，后边改进后续研究了再更新）：

2.3.1 使用2-piece affine gap cost比对

首先了解参考文献中原始算法[Gotoh,O. (1982) An improved algorithm for matching biological sequences.J. Mol. Biol., 162, 705–708.]：
P表示在列水平上，前边分数对当前分数的影响；Q表示在行水平上，前边分数对当前分数的影响。都可以算出具体罚分的，取最小的分数，这里是找出最小的加分，minimap2论文中是找出最大罚分的分数，一个通过加分算的，一个通过减分算的，所以一个是min一个max。

举个栗子，以上边红框为例：

算法如下：

minimap2对参考文献中的算法进行了改进，这里它根据gap的长度对其罚分进行了不同的处理，分为了短gap和长gap对其进行了不同的罚分，其中E就是前面列对当前分数的影响，也就是前面的P；F是前面行对当前分数的影响，也就是前面的Q。通过找到最大的分数得到最优比对结果。
（后面还有对其的进一步改进，后面有空进一步再对其进行解释，还没看:( ）

2.4 剪接序列的对齐

chaining gap cost计算：

这里在上面罚分的基础上，添加了a(j)和d(i)两个受体和供体的罚分，来处理剪接位点。

以经典的剪接位点对其进行不同的罚分。

3 结果

3.1 长基因组读段序列的对齐

3.2 长拼接读段的对齐

3.3 短基因组读段序列的对齐

3.4 长读段装配的对齐

总结

Minimap2是一个多功能的核苷酸序列映射器和成对对准器。

1. 它可用于短读、装配群叠和长噪声基因组和RNA-seq读段，并可作为读段映射器、长读重叠器或全基因组对准器。
2. Minimap2也是精确和高效的，通常在速度和精度方面优于其他领域特定的对齐工具。
3. 现代主流对准器通常使用全文索引来索引引用序列，例如后缀数组或FM-index。这种方法的一个优点是我们可以使用任意长度的精确种子，这有助于增加种子的唯一性并减少不成功的扩展。Minimap2索引使用散列表引用k-mers。这种定长种子在理论上不如变长种子，但实际计算效率更高。当一个查询序列有多个种子命中时，我们可以跳过高度重复的种子，而不会影响最终的准确性。这进一步减轻了对种子唯一性的担忧。同时，在低序位同一性下，长种子也很少见。哈希表是映射长噪声序列的理想数据结构。

[查找的资料]:

Limitations：

• Minimap2 may produce suboptimal alignments through long low-complexity regions where seed positions may be suboptimal. This should not be a big concern because even the optimal alignment may be wrong in such regions.
  Minimap2可以通过种子位置可能是次优的长低复杂度区域产生次优排序。这应该不是一个大问题，因为即使是最优对齐也可能在这些区域出错。
• Minimap2 requires SSE2 instructions on x86 CPUs or NEON on ARM CPUs. It is possible to add non-SIMD support, but it would make minimap2 slower by several times.
  Minimap2需要x86 cpu上的SSE2指令，或者ARM cpu上的NEON指令。可以添加非simd支持，但是这会使minimap2慢几倍
• Minimap2 does not work with a single query or database sequence ~2 billion bases or longer (2,147,483,647 to be exact). The total length of all sequences can well exceed this threshold.
  Minimap2不能处理单个查询或20亿个碱基或更长的数据库序列(准确地说是2,147,483,647)。所有序列的总长度都可能超过这个阈值。
• Minimap2 often misses small exons.（有一篇文章针对这个做了，预发表的Accurate spliced alignment of long RNA sequencing reads（写的条理比较清晰，感觉比minimap强） ）
  Minimap2经常漏掉小的外显子。

流程：

首先，将基因组序列的minimizer存储在哈希表中（minimizer指一段序列内最小哈希值的种子）；
然后，对于每一条待比对序列，找到待比对序列所有的minimizer，通过哈希表找出其在基因组中的位置，并利用chaining算法寻找待比对区域；
最后，将非种子区域用动态规划算法进行比对得到比对结果。

• 搜索minimizer
  minimizer指的是一段序列内最小哈希值的种子，也就是哈希值最小的k-mer。k-mer是长度为k的序列子片段。DNA序列由A、C、G、T四个字符组成，按照计算机编码可以看成一个四进制数。那一个k-mer就可以看做k位的四进制数。比如GCT的哈希值就是2×4的2次方+1×4的1次方+3×4的0次方=39，所以GCT的哈希值就是39。那么可以算出每一个k-mer的哈希值，取w窗口内最小哈希值的k-mer，就是作者定义的minimizer。

  minimap2首先计算基因组序列的minimizer，存储到哈希表中。然后计算待比对序列的minimizer，通过哈希表就可以查找与基因组中一样的minimizer在基因组中的位置。这样每一个minimizer包含三个信息：（1）在基因组中的位置；（2）在待比对序列中的位置；（3）minimizer长度。

• chaining算法：
  chaining就是从上面寻找的一组minimizer集合找出一组共线性的minimizer。chaining方法类似于序列比对过程中的动态规划算法，主要用于找到一组比对区域。共线性就是在基因组中的位置时从左到右排列的。因为相似的序列，肯定包含一些列相同的minimizer。而且序列间越相似，含有相同的minimizer就越多。chaining就是可以找到序列间含有minimizer密度最高的区域，方便后续的比对。

• 种子-扩张阶段：
  通过chaining就找到一组minimizer后，一个minimizer就是一个种子，也是待比对序列和基因组匹配的区域，下一步只需将序列的非种子区域进行比对，与种子区域连接起来，就是最后的序列比对结果。类似于BLAST思想。非种子区域一般比较短，当然是相对整条待比对序列来说的。这样就可以用传统的NW算法或者SW算法进行比对。

[以上仅是自己的理解，欢迎批评指正！]