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

Minimap2：pairwise alignment for nucleotide sequences

前言
1. 引言
2. 方法
- 2.1 计算minimizer
- 2.2 chaining
- - 2.2.1 找到最优chaining scores
  - 2.2.2 回溯
  - 2.2.3 识别primary chains
  - 2.2.4 评估每个序列的散度
  - 2.2.5 均聚物压缩k-mers进行indexing
- 2.3 对齐基因组DNA
- - 2.3.1 使用2-piece affine gap cost比对
- 2.4 剪接序列的对齐
3 结果
- 3.1 长基因组读段序列的对齐
- 3.2 长拼接读段的对齐
- 3.3 短基因组读段序列的对齐
- 3.4 长读段装配的对齐
总结

前言

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

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

1. 引言

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

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

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

2. 方法

根据参考基因组建立minimizer索引表：将基因组序列的minimizer存储在哈希表中（minimizer指一段序列内最小哈希值的种子，因为所研究问题是将RNA-seq比对到参考基因组，因此，是对参考基因组建立索引）；
对于每一条待比对序列（也就是三代测序的reads读段），找到待比对序列所有的minimizer，根据minimizer[x-w+1,x]中x的位置建立anchor表；
根据最优chain分数，将anchor连接起来，从而在参考基因组上进行定位；
再对非种子区域利用动态规划进行扩增。
[注：这里参考基因组建立了一个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是一个多功能的核苷酸序列映射器和成对对准器。

它可用于短读、装配群叠和长噪声基因组和RNA-seq读段，并可作为读段映射器、长读重叠器或全基因组对准器。
Minimap2也是精确和高效的，通常在速度和精度方面优于其他领域特定的对齐工具。
现代主流对准器通常使用全文索引来索引引用序列，例如后缀数组或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算法进行比对。

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

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

SIFT特征点提取及描述论文算法详解
SIFT特征点提取及描述论文算法详解 1. 尺度空间极值检测(Scale-space extrema detection) 1.1 尺度空间和极值 1.2 DoG和LoG的关系 1.3 构建高斯尺度差 ...
联邦学习入门（二）-Practical Secure Aggregation for Privacy-Preserving Machine Learning论文算法详解
本文介绍了一种实用的安全聚合算法(Secure Aggregation Protocol),主要参考了Google的论文Practical Secure Aggregation for Privacy ...
DES（Detection with Enriched Semantics）算法详解
DES算法详解论文背景实验结果算法简介算法结构语义信息增强过程与之前算法的区别算法细节 Semantic enrichment at low level layer Semantic e ...
蚂蚁算法python_Python编程实现蚁群算法详解
简介蚁群算法(ant colony optimization, ACO),又称蚂蚁算法,是一种用来在图中寻找优化路径的机率型算法.它由Marco Dorigo于1992年在他的博士论文中提出,其灵感 ...
PnP算法详解（超详细公式推导）
PnP算法详解 PnP概述 PnP数学模型 PnP求解方法 DLT直接线性变换法 EPnP EPnP的特点步骤理论推倒 1.控制点及齐次重心坐标系 2.控制点的选择 3.计算控制点在相机坐标系下的 ...
人脸识别系列三 | MTCNN算法详解上篇
前言我们前面分享了PCA,Fisher Face,LBPH三种传统的人脸识别算法,Dlib人脸检测算法.今天我们开始分享一下MTCNN算法,这个算法可以将人脸检测和特征点检测结合起来,并且MTCNN ...
Scale-Transferrable Object Detection算法详解（基于多尺寸目标检测）
Scale-Transferrable Object Detection算法详解论文背景算法背景算法简介算法对比算法详解网络结构 DenseNet STM Object Localizat ...
YOLOv4算法详解
YOLOv4: Optimal Speed and Accuracy of Object Detection-论文链接-代码链接目录 1.需求解读 2.YOLOv4算法简介 3.YOLOv4算法详解 ...
目标检测 RCNN算法详解
原文:http://blog.csdn.net/shenxiaolu1984/article/details/51066975 [目标检测]RCNN算法详解 Girshick, Ross, et al ...

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