BWA-MEM算法记录
- BWA-MEM算法总体流程以下：
1 读入 bwt、options、reads；算法
2 利用mem_chain生成chain；数组
3 利用mem_chain_flt过滤掉部分chain；数据结构
4 利用mem_chain2aln生成比对结果元数据。优化
BWA采用seed-and-extend策略。在seed阶段，BWA取read的碱基片断在reference上进行精确匹配，并选择知足必定匹配次数和长度要求read片断做为seed，这个阶段算法的核心是基于FM-index的精确匹配；在extend阶段，BWA利用Smith-Waterman算法将seed在read和reference上向两边延伸比对（容忍gap），进而找到整个read在reference上符合条件的全局匹配。排序
- BWT变换
BWT算法的主要流程是：首先将给定的字符串每次右移一个字符，获得的全部字符串构成一个字符串矩阵，而后对该字符串矩阵中的全部字符串执行字典排序，取该字符串矩阵的最后一列即为BWT变换字符串。下图为对字符串T=“acaacg”进行BWT变换（a）、BWT逆变换（b）、查找字符串（c）的过程。索引
- FM-index数据结构
FM-index是基于BWT的后缀查找算法，时间复杂度为O（m），m为待查找字符串长度。基于FM-index的序列比对算法内存占用较小。内存
下图为字符串T（“acaacg”）构建SA（Suffix Array，后缀数组）和BWT变换的过程。字符串
FM-index主要包含SA数组、OCC数组和BWT索引串三个辅助数据结构。SA数组保存了BW矩阵中后缀字符串在原始串中出现的位置坐标。OCC[c,r]标志BWT（T）中第r行以前字符c的出现次数。it
sp和ep分别表明待匹配序列在SA数组中出现的起始和结束位置（BW矩阵的行号）。io
查找的过程就是将待查找字符串从后向前逆序扩展，利用LF映射的性质进行匹配，延伸至待查找字符串的第一个字符时结束，此时区间值[sp，ep]就是在BW矩阵上精确匹配的区间。区间宽度ep-sp+1表明待查询序列在T中的匹配个数。
LF映射有两个特色：（1）同一行中 F的字母必定是L的后一个字母；（2）相同字母在排序后相对位置不变（例如：F列的第5个a必定是L列的第5个a）
- FM-index存储优化
当reference比较长时，FM-index的内存占用是至关可观的，仅仅OCC数组可能达到几个GB。BWA主要思想是：BWA只存储SA后缀数组和OCC数组的一部分，其余未存储部分则在须要的时候即时计算出来。
BWA的OCC[a,i]只存储i为128的整数倍元素，知足这个条件的i称为一个checkpoint。若是i不是检查点时，首先计算出离i最近的检查点j，而后从OCC数组取出OCC[a,j],最后加上A从B[j]到B[i]出现的次数，结果就是a在BWT变换的位置i以前出现的次数。
SA[i]的空间压缩：当i等于32的整数倍时，后缀数组的元素S[i]会被记录下来，当待查找的序列结束位置不被SA存储时，继续作sp和ep的查找直至找到存储在SA数组中的某一元素，此时用这个元素的值减去该查找步骤进行的步数，便可还原查找串的真实位置。
- 基于SMEM的seed生成
MEM是一段不能继续正向或逆向扩展的精确匹配，SMEM是一种MEM，不被其余MEM所包含。BWA-MEM算法要求对于查询序列的任一个位置，覆盖该位置的最长精确匹配必须是一个SMEM。使用FMD-index和与之相关的正向-逆向扩展算法，能很快地找到全部的SMEM，SMEM查找算法以下图：
BWA-MEM使用传统的seed-and-extend范式，算法中以SMEM做为seed来进行后续的比对。seed是read序列的一部分，用于在比对的初始阶段衡量匹配程序，有的时候，SMEM并不出如今最终的比对结果中。为了减小因为SMEM选择不当致使的误比对，有时候须要从新生成种子（re-seeding）。好比，假设咱们找到一个长度l的SMEM，它在reference中出现了k次。当l太大时，须要从新生成种子，BWA-MEM会找出覆盖原SMEM中间位置碱基的最长精确匹配做为种子，而且要求该种子在reference中至少出现k+1次，新种子的生成算法只须要对SMEM的生成算法进行少量修改。
生成的seed会被后续过程作链化处理。一组距离相近且共线的seed被称为一条链chain，全部的seed链会按照必定的规则作一次筛选。若是一条链的长度小于某一个阈值，或者该链被其余更长的链所包含，则这条链会被抛弃。（该条链的seed也被舍弃？）
算法在seed阶段对于FM-index的访问和计算开销超过总耗时的50%，主要瓶颈是对于FM-index数据结构中OCC矩阵和后缀数组SA的访问。访问OCC数组是为了计算待查找字符串在BW矩阵的匹配区间，访问SA后缀数组则为了定位精确匹配的子序列（seed）在目标序列的起始位置。
- seed扩展
BWA-MEM使用带状动态规划算法对全部seed进行扩展。在对seed进行扩展前，依据链的长度和seed的长度对全部的seed进行排序。若是一条链长度较大，则属于这条链的全部seed获得的排序就比较靠前，基于链长的排序结束后，再按照seed长度从长到端对链中全部的seed作一次排序。每个seed按照排序的前后进行扩展，首先判断它是否被以前的比对结果所包含，若是没有被包含，则使用带状动态规划算法对它进行扩展。BWA-MEM使用的带状动态规划算法经过对Smith-Waterman算法进行SSE2指令加速。
Banded DP的扩展过程当中，有可能会遇到这样的一段序列：它的中间部分比对结果不好，而它的两侧部分的匹配质量很高。对于这样的区域，BWA-MEM使用了Z-dropoff的启发式思想，当比对的分数从最大值降低到低于某个阈值时，将会中止序列比对过程。
BWA-MEM算法会记录比对到query末端的最大得分，当这个分数和局部对比的分数的差小于某一个阈值时，即便局部比对的分数较高，BWA-MEM也会将它丢弃而采用比对到末端的得分。BWA-MEM的总体算法流程以下图。（？？？）
- OCC数组
OCC数组是按照必定的间隔存储的，在BWA-MEM算法中，OCC数组每隔128bp存储一个元素。在经典的FM-index结构中，OCC数组和BWT数组是分别存储的，而BWA-MEM的FMD-index中，这两个数组被紧凑地合并到了一个数组中，OCC和BWT组成的新数组的基本单元（一个cell）是由OCC数组和reference的BWT变换数组拼成的64个字节，其中前32个字节分别存储在BWT数组的位置i*128以前ATCG四种碱基出现的总次数，即四种碱基对应的OCC数组元素，后32个字节用来存储BWT数组在位置i*128以后128bp的碱基。