比对算法总结(二)——基于BWT索引结构的比对算法-Bowite1

Posted 2023-03-30

参考技术A

这是美国马里兰大学计算机研究所、生物信息学和计算生物学中心于2009年发表在《Genome Biology》杂志的一篇经典文章，至此以后依赖于BWT索引的比对算法成为主流。 Bowite 是一款超快速、内存占用低的短序列比对软件，适用于将短reads比对至大型参考基因组。采用Burrows-Wheeler 算法建立索引的Bowite软件可以在1 CPU时内，将2000万条reads 比对至人参考基因组，且内存只占有1.3Gb。于此同时Bowite 采用了新的quality-aware backtracking（质量回溯）算法，比对过程允许错配。

在此之前都是采用对reads (SHRiMP, Maq, RMAP,ZOOM) 或者参考基因组 (SOAP)构建哈希表的算法进行序列比对，该算法已在上篇文章中进行了介绍 https://www.jianshu.com/p/f5ccff73b181 。
Bowite 采用了一种完全新的索引构建策略，适用于哺乳动物重测序。根据千人基因组计划数据，Bowite 在35bp PE 序列上的比对速度要比Maq 软件快35 倍，比SOAP软件快300倍。Bowite 采用 Burrows-Wheeler 算法对 full-text minute-space (FM) 构建索引，人参考基因组占用的内存为1.3 GB。
为了追求速度，Bowite 针对哺乳动物重测序项目进行了很多合理的折中。例如，如果一条reads有多条最优匹配，Bowite 只会输出一条最优匹配。当输出的最优匹配也不是完全匹配时，Bowite并不能保证在所有情况下都能输出最高质量的匹配。在设定了较高的匹配阈值时，一小部分含有多个错配的reads可能会比对失败。在默认参数条件下，Bowite 的灵敏度与SOAP 相当，略低于Maq。可以在命令行手动改变参数，在牺牲更多时间的情况下，增加灵敏度，给出reads所有可能的比对结果。目前Bowite 比对的reads长度范围为4bp - 1024bp。

Bowite 对参考基因组建立索引的方法是 Burrows-Wheeler transform (BWT) 和 FM index。Bowite 建立的人类基因组索引在硬盘上的大小为2.2GB，在比对时的内存为1.3GB。FM index 常用的精确查找方法为 Ferragina 和 Manzini 算法。Bowite 没有完全使用该算法，因为该算法不允许错配，不能比对含有测序错误和变异的reads。针对这种情况，Bowite引入了新的扩展算法：quality-aware backtracking 算法，允许错配并支持高质量比对；double indexing 策略，避免过度回溯；Bowite比对策略与Maq软件相似，允许小部分的高质量reads 含有错配，并且对所有的错配位点的质量值设置了上限阈值。

BWT 转换是字符串的可逆性排列，它最早应用于文本数据的压缩，依赖BWT建立的索引，可以在较低内存下，实现大型文本的有效搜索。它被在生物信息学中有广泛的应用，包括重复区域计数、全基因组比对、微阵列探针设计、Smith-Waterman 比对到人参考基因组。Burrows-Wheeler transform (BWT) 的转换步骤如图1所示：

1、轮转排序。如图1a 所示，（1）将字符$ 添加到文本 T （acaacg）的末尾，但需注意其中字符$ 并未实际添加到文本 T 中，且其在字母表中逻辑顺序小于 T 中所有出现过的字符。（2） 然后将当前字符串的第一个字符移到最后一位，形成一个新的字符串，再将新的字符串的第一位移到最后一位形成另一个新的字符串，就这样不断循环这个过程，直到字符串循环完毕（即$处于第一位）,这样就形成了一个基于原字符串的字符矩阵M（这一步原图1a 进行了省略，见下方小图）。（3） 然后对矩阵M的各行字符按照字典先后顺序排序，获得排序后的字符矩阵 BWM（T），矩阵的最后一列定义为 BWT（T）。 前期经过一个小复杂的过程获得了BWT（T）列，那这一列到底有什么用呢？其实BWT（T）列通过简单的算法就可以推算出原始文本T的所有信息。而经过转换之后的BWT（T）列大量重复字符是靠近的，只储存该列信息，可以大大提高字符压缩比例。

2、LF-Mapping。图1a 转换矩阵 BWM（T）含有一种 \'last first (LF) mapping\' 的特性，即最后一列L中出现某字符出现的顺序与第一列F某字符出现的次序时一致的。根据Supplementary1 图中算法1 STEPLEFT 和 算法2 UNPERMUTE 就可以推算出BWT（T）到 T 的过程， 图1 b记录了整个推算过程。 详细推算过程可参考这个博客介绍： https://blog.csdn.net/stormlovetao/article/details/7048481 。

3、reads精确匹配。使用BWT算法的最终目的是要将短reads比对到参考基因组上，确定短reads在参考基因组上的具体位置。转换后的BWT（T）序列，可以利用Supplementary1 图中算法3 EXACTMATCH 实现reads的精确匹配。图1c 列出了 字符串 aac 比对至acaacg 的过程 。 详细推算过程可参考这篇介绍： https://zhuanlan.zhihu.com/p/158901556 。

上述的BWT转换只能用于精确的匹配，但是测序reads是含有测序错误和突变的，精确匹配并不适用。这里应用了 backtracking 搜索的算法，用于允许错配快速比对 。含有错配的reads只是一小部分。测序reads的每个碱基都含有唯一的测序量值，测序质量值越该位点是测序错误的可能越大，只有当一条read 的所有错配的测序质量值总和小于一定阈值时可以允许错误匹配。
图2显示了精确匹配和非精确匹配的过程，backtracking 搜索过程类似于 EXACTMATCH ，首先计算连续较长的后缀矩阵。如果矩阵中没有搜索到相应的reads，则算法会选择一个已经匹配的查询位置，替换一个不同碱基，再次进行匹配。EXACTMATCH搜索从被替换位置之后开始，这样就可以比对就可以允许一定的错配。backtracking 过程发生在堆栈结构的上下文中，当有替换产生时，堆栈的结构会增长；当所有结果都不匹配时，堆栈结构会收缩。
Bowite 软件的搜索算法是比较贪婪的，Bowite软件会报出遇到的第一个有效比对，并不一定是在错配数目和变异质量上的“最佳比对”。没有查询最优比对的原因是寻找“最佳比对”会比现有的模型慢2-3倍。而在重测序项目上，速度是更重要的因素。Bowite 也设置了可以输出多个比对位置（-k）和所有比对位置（-a）的参数，添加这些参数后，比对速度会显著变慢。

目前的比对软件会有过度回溯的情况，在reads的3‘端花费大量无用时间去回溯。Bowite利用‘double indexing’技术减少了过度回溯的发生。简单来说就是对正向参考基因组进行BWT转换，称为 ‘Forward index’，同时对反向（注意不是互补配对序列,是反向序列）参考基因组也进行BWT转换，称为‘Mirror index’。 当只允许一个错配时，比对根据reads是前半段出现错配，还是后半段出现错配会有两种情况：（1）Phase1 将Forward index 加载入内存，不允许查询reads右半段出现错配；（2）Phase2 将Mirror index 加载如内存，不允许查询序列的反向reads右半段（原查询序列的左半端） 出现错配。这样可以避免过度回溯，提高比比对的灵敏度。 但是，如果比对软件允许一个reads有多个错配时，仍然会有过度回溯的现象发生，为了减少过度回溯现象的发生，这里将回溯的上限进行了限定（默认值为：125次）。

Bowite 允许使用者在高质量reads的末端（默认是28bp）设置错配数目（默认的错配数目是2）。高质量reads末端的28bp序列被称为 \'种子\' 序列。这个‘种子’序列又可分为两等份：14bp的高质量末端称为 ‘hi-half’（通常位于5‘端）,14bp的低质量末端称为‘lo-half’。 如果种子序列只允许2bp 的错配，比对会出现4 种情况：（1）种子序列中没有错配(case1)；（2）hi-half区域没有错配，lo-half区域有一个或两个错配（case2）;（3）lo-half区域没有错配，hi-half区域有一个或两个错配（case3）;（4）lo-half区域有一个错配，hi-half区域有一个错配（case4）；
在所有情况下，reads的非种子部分允许任意数目的错配。如图3所示，Bowite 算法会根据上面4 种情况交替变化‘Forward index’和‘Mirror index’比对策略，主要会有三种比对策略。

Bowite 建立一次参考基因组索引后，后续的比对可反复使用该索引。表1和表2列出了在默认参数条件下，Bowite、SOAP、Maq软件性能的比较。在reads比对率相近的条件下，Bowite软件的比对速度速度相对于SOAP、Maq软件有较大的提升。

1、将reads 比对至人参考基因组上，Bowite相对于SOAP和Maq软件有较大的优势。它运行的内存非常小（1.2GB），在相同灵敏度下，速度有了较大的提升。
2、Bowite 软件建立一次参考基因组索引后，后续的比对可反复使用该索引。
3、Bowite 速度快、内存占用小、灵敏度高主要是因为使用了BWT算法构建索引、利用回溯算法允许错配、采用Double index策略避免过度回溯。
4、Bowite 软件目前并不支持插入、缺失比对，这个是今后需要努力的方向。

[1] Langmead B . Ultrafast and memory-efficient alignment of short DNA sequences to the human genome[J]. Genome biology, 2009, 10(3):R25.
[2] BWT 推算过程参考博客 https://blog.csdn.net/stormlovetao/article/details/7048481
[3] FM index 精确查匹配过程参考文章 https://zhuanlan.zhihu.com/p/158901556

基于BWT算法的比对软件原理解析（BWA & Bowtie & Bowtie2）

参考技术A

参考：
踏踏实实做技术：BWA，Bowtie，Bowtie2的比对算法推导

remove multiple mapping reads的方法

CHIP-seq: Bowtie2、BWA用的比较多
RNA-seq: Tophat、Bsmap
甲基化：BS-seeker

global---NW
local--SW

好处是能够穷举出所有的比对情况，所以可以选择全局最优的结果；最大的缺点是比对的非常慢。

BWT（Burrows-Wheeler Transform ）

第一步，在raw seq中加$符号，并平移，形成一个 raw matrix

第二步，根据Raw Matrix的首字母进行排序，得到转换矩阵Matrix’，默认$符号排在第一位，

所以最后只用保存L列和每个字母的相对位置就可以了，根据L列和每个字母的相对位置可以干两件事情：

例如：第一个是L- 对应F- 的前一个是G，L-G对应F-G；F-G的前一个是L-C，依次类推，得到原来的ref：ACAACG$

14bp（high quality）---14bp（low quality of high quality）--8bp（real low quality）
分成三断seed，seed1+seed2比对总共的mismatch <= 2，则继续8bp的比对；如果 > 2 直接放弃后面的比对；

第一步，选择seed区域；
20里面选18---
（18+2）+（18+2）+（18+2）+...+(18+2)
保证一个fragment是20，seed 是18bp
或者，10里面选16--
fragment = 16，overlap = 6，

那么根据BWT算法，就把拆分的seed mapping到基因组的大概位置；
然后把基因组可能mapping上的那段区域挑出来，和query seq做比对（用NW或者SW算法），因为query seq NW和SW允许gap open