seed alignment 算法（BWT）

人类的参考基因组大约含有31.6亿个碱基对，100万条150bp读长的read假若一一对应得通过Needleman-Wunsch算法（https://www.jianshu.com/p/85465f32273a）做双序列比对，时间真让人等不起。因此，大多数mapping软件将“序列回帖”分成 seed alignment 和 seed extension两个步骤。
seed alignment：首先通过截取read中的短序列片段（称为 seed）与参考基因组比对，来找到reference的index；
seed extension：然后通过reference的index将附近的序列与seed对应read做双序列比对 （Needleman-Wunsch算法或Smith-Waterman算法）

这篇文章就记录一下seed alignment 的BWT算法，不过，BWT算法原本用于数据压缩，而它的 压缩 和 解压缩 的过程也可以直接类比到，参考基因组做索引 和 seed alignment 的双序列比对。

1. BWT算法做数据压缩的原理

1.1 序列ababc数据压缩

图1.BWT算法压缩原理

以下步骤与图1一一对应：
举例：压缩字符串ababc
输入字符串 ababc
第一步，添加标记ababc$
第二步，ababc$“循环转移”（序列最后一个字母“依次”移动到最前端）
第三步，将“循环转移”获得的矩阵按照 第一列首字母 排序获得M数组
第四步，取出M数组的第一列为 F列；M数组的最后一列为 L列
数据压缩：做到这一步之后，便可以直接将 L列 以c,$,b,2a,b的形式存储，实现了字符串ababc的数据压缩（不过这个举例里压缩率并不高😓）

1.2 数据解压缩

数据解压缩，也就是从L列内容还原原始字符串。其中势必用到了M数组中的两个特殊列 F列和L列，还有他们之间的相互关系。

BWT算法的三项原则：

1. 同一行中，“L列字母”为“F列字母”的前一个字母
2. L列字母排序就可以得到F列（因为F列和L列其实元素是相同的）
3. 单字母的相对位置不变（F列的第n个a对应L列的第n个a，这个结论可以自己推导）

如何通过F列和L列还原原始序列？也就是解压缩过程：

通过L列和F列，反推导原始序列
$为序列的结尾标志，因此从$开始反向推导原始序列（“L列字母”为“F列字母”的前一个字母；F列和L列单字母的相对位置不变）

• 从F列的$开始
• 根据第一行，$的前一个字母为 c
• L列第一行的c对应F列的第六行的c（均为两列中的第1个c）
• 根据第六行，c的前一个字母为 b
• L列第六行的b对应F列的第五行的b（均为两列中的第2个b）
• 根据第五行，b的前一个字母为 a
• L列第五行的a对应F列的第三行的a（均为两列中的第2个a）
• 根据第三行，a的前一个字母为 b
• L列第三行的b对应F列的第四行b（均为两列中的第1个b）
• 根据第四行，b的前一个字母为 a
• L列第四行的a对应F列的第二行a（均为两列中的第1个a）
• 根据第二行，a的前一个字母为 $（已经找回到起点 $，说明序列已经完全恢复）

以上步骤中，黑粗体描述的字母 从下向上 排列为：ababc，即达到了恢复原始序列的目的。

BWT用于数据压缩，是因为L列（c,$,b,a,a,b）可以直接储存为（c,$,b,2a,b）。这种压缩重复碱基的形式，对于存在“高重复区域”的基因组序列而言，可以极大的减少存储空间。同时借助F列做位置参考，很容易能找回原始序列。

2. BWT算法做序列比对

我们如何用BWT的算法做碱基序列比对？实际上，以上提到 数据压缩和解压缩 的过程就是我们做序列比对的过程。
1）我们建立参考基因组的索引，其实便是建立refercen序列的L列和它相对位置的index（体现在👆便是ababc获得L列的过程，也就是 数据压缩 的过程）；
2）我们将测序得到的reads与参考基因组比对，其实便是查找reads对应参考基因组的位置，并观察reads序列是否可以还原出对应位置的碱基序列（体现在👆便是由L列排序获得F列，然后以F列配合做指引，从最后一个字母出发做 数据解压缩 ）

举例：abab是否为ababc的子序列？我们看BWT算法是如何判断的

abab是否为ababc的子序列
比对的过程是从后向前的"倒序比对"，因此abab比对顺序应该为baba，从F列对应字母开始向前追溯。因为F列有2个b，因此要做两次比对：
第一次：

• 从F列的第1个 b 开始
• 根据第四行，b的前一个字母为 a
• L列第四行的a对应F列的第二行的a
• 根据第二行，a的前一个字母为 $
  直接终止比对过程，比对到的序列为ab

第二次：

• 从F列的第2个 b 开始
• 根据第五行，b的前一个字母为 a
• L列第五行的a对应F列的第三行的a
• 根据第三行，a的前一个字母为 b
• L列第三行的b对应F列的第四行的b
• 根据第四行，b的前一个字母为 a
  至此已经比对到序列 abab
  因此软件判断abab为ababc的子序列

PS：bowtie 里的FM-index 的实现方式：
在 BWT 算法的基础上增加以下两个问题的解决方法

⚠️ 1. 计算机如何判断参考序列（即L列）中 碱基的相对位置？
例如：当前seed比对到参考序列的 A 在参考序列（即L列） “所有A 碱基”中的排序号码是多少？
bowtie的做法：每隔128行设置一个check point（A:233;G:231;C:256;T:271），通过距离该碱基最近的check point获得碱基的相对位置。

⚠️ 2. 计算机如何判断reads在参考基因组中的位置？
最简单粗暴的方式是：通过将后缀数组载入内存，添加每个碱基的位置信息。但是人类基因组的后缀数组约为12G，完全载入特别耗内存。
bowtie的做法：每32行设置一个sample位置值，因此要确定位置信息，只需要回读到sample位置区，即可推导read比对到参考基因组的位置信息。（如此便将后缀数组的内存降低为1/32）
参考链接：https://blog.csdn.net/stormlovetao/article/details/7048481

实际的比对过程中，测序得到的reads都被分割成几十bp的片段，选取其中的部分质量较好的序列作为seed序列与参考基因组比对（循环如上的 解压缩 过程），找到reads在基因组上大概的位置（比对上的位置可能会很多，会综合很多因素：insertion、deletion、mismatch、reads quality等等，为每一个位置打分，最终取得分最高的位置）。
确定位置之后，取出参考基因组对应位置附近的序列，和reads做双序列比对。