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
前言
最近研究将三代测序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,从而定位到参考基因组上。以上很多名词后面都会进行一一解释和计算]
[图片仅供参考,感觉不是很准确]
![](https://i-blog.csdnimg.cn/blog_migrate/a7cd07fa1f6e3f70c6d008a28cc6182f.png#pic_center)
主要流程如下:
![](https://i-blog.csdnimg.cn/blog_migrate/798bc92ad3f9decc8bdaca68f2386129.png#pic_center)
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算法进行比对。
[以上仅是自己的理解,欢迎批评指正!]