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

工欲善其事 必先利其器

前言

关于比对软件BWA前面我们介绍了

• 序列比对之BWA
• 比对软件BWA及其算法（上）

今天再来细说一下BWA-MEM

一、BWA-MEM2介绍与下载使用

BWA-MEM是李恒大神于2010在bioinformatics发布的一款比对软件

Fast and accurate short read alignment with Burrows–Wheeler transform: https://doi.org/10.1093/bioinformatics/btp324。

BWA-MEM2是英特尔并行计算实验室和李恒大神于2019年发布的，运用C++并行计算技术和更高效的处理器指令集加速BWA-MEM算法的比对软件。

Efficient Architecture-Aware Acceleration of BWA-MEM for Multicore Systems: https://ieeexplore.ieee.org/document/8820962

BWA-MEM2源代码和软件下载：bwa-mem2/bwa-mem2: The next version of bwa-mem (https://github.com/bwa-mem2/bwa-mem2)

#方法一conda安装：
conda install bwa-mem2

#方法二下载源代码或者压缩包安装，进入文件夹下进行编译
make #make编译要C语言编译器，sudo yum install gcc。若报错也可能是下的源代码中safestringlib是空的，确保自己下的bwa-mem2的bwa-mem2-2.2.1/ext/safestringlib文件夹不为空

二、BWA-MEM使用

在进行比对之前我们需要构建参考基因组的索引，这个过程类似于建立字典的目录，以方便快速查询字典一样。构建索引是比较占用内存的。构建人参考基因组索引建议服务器（或自己的工作站）的内存至少有128G。

#查看使用方法：
bwa-mem2 index
#bwa-mem2 index [-p prefix] <in.fasta>
#对参考基因组test.fa构建索引
bwa-mem2 index test.fa #当-p参数缺省时，会默认将参考基因组文件名作为前缀生成索引文件

-p参数是生成的索引文件的前缀，in.fasta参数是fasta格式（可以gz压缩）的参考基因组。

他会产生以下后缀的索引文件（做个了解，其实会用就行）：

.0123（二进制文件，0123对应ACGT的参考基因组）

.amb（参考基因组上连续的N（holes）的位置）

.ann（参考基因组名称，长度，起始位置等信息）

.bwt.2bit.64（二进制文件，bwt和后缀数组等信息）

.pac（二进制文件，两位（一个字节中存四个碱基）代表一个碱基的压缩版参考基因组）

进行比对bwa-mem2 mem

#查看比对参数，使用方法
bwa-mem2 mem
#Usage: bwa-mem2 mem [options] <idxbase> <in1.fq> [in2.fq]

除这些最基本的参数外，还可以使用其他参数自定义程序使用线程数，种子序列长度，比对打分罚分等等

这里还没来得及整理，以后会根据源代码仔细讲讲每个参数用法。这里我们使用最基本的参数和-o参数指定输出文件（若不使用-o参数会把sam文件打印到屏幕上）就好。

bwa-mem2 mem [options] \<idxbase\> <in1.fq> [in2.fq]中：options是上图展示的可选参数，不选时都使用默认值；idxbase则是我们构建索引时-p参数生成的前缀名；in1.fq是单端测序的fastq或者双端测序的fastq1；in2.fq是可选参数，为双端测序的fastq2

#我们使用上文对ref.fa构建的索引，比对双端测序数据test.fq1 test.fq2
bwa-mem2 mem ref.fa test.fq1 test.fq2 -o test.sam

这样我们就完成比对得到sam格式比对结果了。

三、BWA-MEM算法简要介绍

我们使用大写字母X,Y等来表示DNA序列，字母表∑=A,C,G,T上的字符表示四种可能的核苷酸（也可称为碱基）。我们用表示序列|x|的长度。设X[i]表示序列X上位置i处的碱基，X[i,j]（i≤j）表示X的片段X[i]X[i+1]X[i+2]。

设Q为一个短读段查询(query)序列，R为一个长的参考(reference)基因组序列。序列比对问题是找出Q在R上的最匹配结果。对于短读段序列，Q的长度范围取决于测序平台，通常在50到250bp之间。|R|可以从几十万（对于细菌基因组）到几十亿bp（对于植物基因组）不等。人类的参考基因组长度为30亿个碱基对。大多数流行的序列比对工具，包括BWA-MEM，使用经典的种子序列定位及延伸算法(seed and extend)。在播种阶段，找到读段的短子字符串（称为种子序列）在参考序列中的精确比对，允许比对中有零或非常少量的差异。这给出了整个读段可能比对到的位置。在延伸阶段，延伸种子序列的两侧直至覆盖整个读段，通常使用基于动态规划的算法如Smith-Waterman算法(Smith and Waterman 1981)，计算每个比对位置的得分并报告最佳比对结果。

3.1 FM索引构建

BWA-MEM使用参考基因组序列的FM索引来生成种子序列。参考基因中的FM索引基于其Burrows-Wheeler变换(BWT, Burrows-Wheeler transform)和后缀数组(SA, suffix array)。图1展示了如何构建示例序列R的BWT、后缀数组(SA, suffix array)。首先，将R的末尾附加上结束字符$，我们认定它在∑中的字母表顺序小于所有字符。随后，获取 的所有旋转(Rotations)序列。对这些旋转序列按字母表进行排序得到BW矩阵。BWT(S)是该矩阵的最后一列。后缀数组(S)存储这些旋转的第一个碱基在R 中的原始位置，即R的后缀的排序顺序。

图1

图1（Fig.1）构建参考基因组的后缀数组和BWT：我们以序列R作为示例参考基因组，图中左侧矩阵是由R经旋转(Rotation)获得，称为旋转矩阵，标红的碱基是参考基因组的第一个碱基。将该矩阵按照字母表顺序Σ排序，得到图中右侧矩阵，称为BW矩阵。BW矩阵第i列在旋转矩阵的原先位置为S[i]，S称为后缀数组(Suffix array)，右边矩阵的每行首位碱基组成F列(First Column)，矩阵的每行末位碱基组成L列(Last Column)。F列是每种碱基按字母表顺序重复其在参考基因组中出现的次数，L列即为BWT字符串(Burrows-Wheeler transform)。

查询读段的所有精确比对都是BW矩阵中旋转序列的前子字符串。因为BW矩阵像字典的索引一样，按字母表顺序排序，所以这些比对会处在BW矩阵的连续行中。因此，查询读段的所有比对可以表示为BW矩阵数行的范围。这个范围被称为查询读段的后缀数组区间(Suffix Array interval)。在图1中，序列“AT”比对到了SA区间[1,2]。为了更快速地检索BWT，我们使用了FM索引(Ferragina-Manzini Index)(Ferragina and Manzini 2001)，如图3、图4。它由D和O矩阵组成。D[x]是在R[O,|R|-1]（不包括$）中字典顺序小于x∈∑的碱基的数量，而O[x,i]是B[0,i]中x的出现次数。

图2

图2（Fig.2）对查询序列的精确检索：在获得F列和L列之后，我们通过LF比对回溯查询序列，这一点之后会在图5中详细解释，回溯比对得到的结果在BW矩阵上时一个区间，称为后缀数组区间(SA, Suffix array interval)，后缀数组区间中的每一个S值都对应到参考基因组上的一个位置。图中查询序列AT比对到SA区间[1,2]，区间中的S值5和0表明在参考基因组R的第5和0个位置开始为AT。

图3

图3.（Fig.3）O矩阵：O[x,i]为碱基x在B[0,i]中累计的个数，表格的列为BWT上的第几个碱基位置，行为四种碱基。O矩阵的作用是压缩L列，通过局部载入后缀数组S和压缩后的L列，实现在内存中对整个后缀数组S进行动态计算，极大的减少了后缀数组S和BWT的内存占用。

图4

图4.（Fig.4）D矩阵：D[x]为在R上按字母表顺序∑小于碱基x的碱基个数（不包含$）。D矩阵的作用是压缩F列，减少其内存占用。

3.2 BWA-MEM算法

BWA-MEM算法是BWA-MEM2软件包执行比对的核心部分，用于将测序得到的读段序列比对到参考基因组上。其主要步骤如下：

3.2.1 SMEM

寻找超精确比对(SMEM, super maximal exact matches)：BWA-MEM首先寻找读段与参考基因组之间的超精确比对。最大精确比对(MEM, maximal exact matches)是读段的子字符串在参考基因组上的精确比对，且不能在任何方向上进一步延伸。超精确比对是查询读段每个位置中覆盖该位置的最长精确匹配。BWA-MEM使用FM索引进行LF比对(Last-to-First Mapping)来高效地检索SMEM，获得SMEM的后缀数组区间，LF比对详细原理如图5。假设我们的查询读段为融合转录本，其靠近5’端的序列TT来自参考基因组染色体1上的某基因，靠近3’端的序列GAT来自参考基因组染色体2上另一基因。因为LF比对是自后向前回溯的，所以我们首先从查询序列3’端的T开始，根据之前旋转(Rotation)的规则，同一行中L列的碱基实际上是F列中的前一个碱基，所以比对从F列中的三个T起始，这三行中有两行F列碱基为T且L列碱基为A，因此我们启动红色和绿色两条比对路径。BW矩阵还有一个性质，即F列中的出现的第k个字符x和L列出现的第k个字符x在原字符串R上是同一个字符。具体来说，红色比对路径L列中的A1 会对应到F列中的A1，我们用红色箭头和下角标标出；而绿色比对路径L列中的A2对应到F列中的A2，使用绿色箭头和下角标标出。随后我们沿红色比对路径，F列中的A1其前一个碱基为G2，比对成功，继续延伸；而绿色比对路径在F列中的A2其前一个碱基为$，不为G且比对到参考基因组起始位置，因此终止绿色比对路径。红色比对路径F列中的G1在参考基因组中的前一个碱基为C1 ，而查询序列为T，不匹配所以红色比对路径也终止。至此红色比对路径产生了一个SMEM（在实际情况中要求SMEM至少大于19bp）因此我们可以接着启动蓝色比对路径，蓝色路径会沿红色路径比对方向重新开始一个比对，该方向比对完之后为寻找SMEM，从示意图虚线位置逆向往回比对，两个比对方向都无法延伸时得到第二个SMEM。若此时查询读段还有部分长度未进行LF比对则会循环蓝色路径的比对步骤直至该读段所有长度都进行了LF比对。示例中产生的两个SMEM其S值如图左下的BW矩阵。

图5（Fig. 5）LF比对回溯查询读段

3.2.2 SAL

执行后缀数组查找(SAL, Suffix Array Lookup)：执行后缀数组查找以获取与前一步中获得的后缀区间在参考基因组中的坐标。由于将后缀数组全部载入占用内存空间过大，所以BWA-MEM通过局部载入后缀数组和F、L列（即经过压缩后的FM索引），通过已经读入内存的后缀数组对剩余未读入的后缀数组进行LF比对回溯(backtrace)，进行动态计算，以实现内存友好的后缀数组查找。

3.2.3 chain

链形成(chain)：BWA-MEM将共线且彼此接近的的种子序列（即筛选过的SMEM）进行链接，如图6。这步链接过程有助于过滤掉假阳性比对的种子序列，并提高了比对的效率。

图6 BWA-MEM算法：步骤1、2为SMEM和SAL示意；步骤3为CHAIN示意；步骤4为BSW示意

3.2.4 BSW

带状Smith-Waterman算法延伸(BSW，banded Smith-Waterman)：使用基于动态规划(DP, dynamic programming)的带状Smith-Waterman比对算法来延伸种子序列(SMEM)。BSW算法仅计算动态规划矩阵的对角线带。这种带状比对方法相比于完整的Smith-Waterman显著提高了比对速度，用于将读段精确地比对到参考基因组中。对于空位罚分，BWA-MEM使用了一种启发式的动态罚分算法，在避免过度延伸的同时也赋予比对程序跨越较差区域的能力。并且当读段末尾出现突变而遭到罚分，而延伸至完整末尾的得分和最优局部比对得分相差小于算法给定阈值时，会拒绝最优比对结果，选择延伸至完整末尾的结果(end-to-end)，保证了末尾出现突变的读段可以正常比对至参考基因组上。

3.2.5 SAM

SAM格式输出：最后，以SAM(Sequence Alignment/Map)格式输出比对结果，该格式是储存序列比对结果的通用格式。