比对软件BWA及其算法（上）

一、BWA与BWT

BWA基础用法参见：序列比对之BWA 官网：Burrows-Wheeler Aligner (https://bio-bwa.sourceforge.net/)

BWA是一个用于将与参考基因组相异性较低的序列比对到较大参考基因组（如人类基因组）的软件。它包含三个算法: BWA-backtrack, BWA-SW and BWA-MEM。

BWA-backtrack 算法是为了illumina测序100bp长的read设计的，其余两个算法用于比对在70到1Mbp之间的较长的序列。BWA-MEM 和 BWA-SW具有相同的特性，如支持长读长和切割比对，但是BWA-MEM这一最新的算法，通常推荐当reads质量较高时使用，因为它更快更准确。且对于70-100bp的Illumina reads，BWA-MEM相比较BWA-backtrack具有更好的性能。

BWA软件在压缩参考基因组，构建参考基因组的索引，以及比对过程中使用BWT算法。BWT算法是M. Burrows和D.J. Wheeler最开始提出对较大字符串文本进行压缩的算法。其部分特性特别适用于我们进行序列的比对。

二、BWT算法

我们以文献中的字符串googol 为例， 代表结束的字符，在字符串中有且仅有一个，且在字母表顺序中排第一位，例如在26字母表中

首先我们要生成左边形式的矩阵，他是将上一行的字符串的第一个字符放到最后一位形成的。随后我们将每一行新的字符串从前到后按字母表顺序排列，生成右边的矩阵，称为Suffix array矩阵，矩阵最后一列 looogg 称为Burrows-Wheeler Transform string （BWT string）。可以看到他将部分o和g重排到一起了，所以我们也可以将这个字符串记为lo2o2g。在这个短字符串的例子中可能无法体现其压缩效率，但是当我们对长字符串如参考基因组处理时，BWT算法可以有效的压缩文本。

BWT算法还有一些特性，我们将SA矩阵的第一列称为F列，最后一列（BWT string）称为L列，明显F列和L列中各字母数量相同，且在原字符串中的顺序相同，如下图所示。

BWT算法是可逆的，即我们知道BWT string和SA矩阵中index为0的字符串，即上图左边矩阵的第一条字符串（我们的原始输入），我们就可以进行backtrace。

在我们进行比对过程中，我们利用SA矩阵将BWT矩阵的string按字母表字典中顺序放在一起的特性，可以像检索字典一样实现快速的比对。例如我们要将序列go比对到参考基因组googol上，那么SA矩阵中所有首字母为g，第二个字母为o的行就是我们比对到的结果。图中红框代表我们比对到的结果，称为go在SA矩阵中的interval，还记得前面3,0的数字代表：图1中左边的BWT矩阵未经过字母表顺序排序时的顺序，它也说明go序列比对到了googol参考基因组的第0个和第三个字符开始的位置，即：googol，googol。

因此我们知道，如果我们考虑无容错的比对时，我们querying的序列必定会比对到参考基因组的一个SA interval中。但是当我们考虑有容错的情况，比如我们允许gg，gl，oo，lo，甚至ao，gb时（这些都是存在一个字符的mismatch），比对上的结果就有可能出现在多个SA interval中，如下图所示。两个红框分别代表了go和oo比对到参考基因组googol上的SA inteval，因为我们允许一个mismatch，所以我们认为oo也是go比对上的结果。