Bowtie 里的FM-index 简介

FM-index是bowtie里的一种索引方式，而bowtie是一种mapping软件，mapping是什么？就是将一个短read 映射到reference genome 中。Read是什么？Reference genome是什么？这个不知道的话，建议就不要继续往下看了。。。。。。

言归正传，FM-index，包括三部分，①BWT（T），②checkpoint data，③一个简化了的SA[]数组。当然在bowtie里面还有其他的一些结构，如lookup，就不一一介绍了。

1.首先说一下BWT（T）

BWT是The Burrows-Wheeler Transform的缩写，是对参考基因组的变换方式，也就是对参考基因组进行了一次有规律的重新排序，变换的目的就是为了方便后续进行查找。

网上对于BWT转化得素材较多，我不赘述，下面通过一个例子简单介绍一下BWT（T）。PS.之所以这么写BWT（T）是因为本身BWT（）是个函数，参数T就是一个字符串，当然在genome mapping里T就是指的参考基因组了。

先从形式上看，左面的acaacg这个串就是我们的参数T，T = acaacg；右面是转换了的T，BWT（T）= gc$aaac。美元符‘$’是我们人为加进去的，定义$的字典序小于其他任何字符。

BWT（T）也叫做轮转排序，通过中间的那一步，我们很容易发现，原串acaacg$通过轮转，形成一个字符串矩阵Burrows-Wheeler Matrices，最后一列就是BWT（T）。记住，BWT（T）就是一个与原串T等长的字符串，只是经过了轮转排序而已。BWT（T）有一个性质是非常重要的，也是理解LF-mapping，FM-index的关键：在Burrows-Wheeler Matrices中，最后一列第i个出现的字符x，与第一列中第i个出现的x是同一个。比如，最后一列中第2个a，与第一列中出现的第二个a是同一个a。细心的同学会发现，第一列其实就是T的一个字典排序。

那这个BWT（T）有什么用呢？

BWT（T）用处之一：首先我们通过BWT（T）可以还原出T，原理就是Last First (LF)Mapping，用到的算法叫做UNPERMUTE algorithm。

先介绍LF-Mapping，其实就是怎么将最后一列的某个字符x映射到第一列中，看这个x在第一列中哪一行。很简单，只需要知道，字典序排在x之前的所有字符的个数，和在最后一列中这个x在所有x中是第几个。因此我们定义了两个参数C[c],和Occ[c,r]。c代表某字符，用上面的例子说，就是acg中的某一个(相当于我之前提到的x，例子中都用c表示字符，我也用c吧。。。)，r代表矩阵的哪一行。C[c]的定义是，字典序小于字符c的所有字符个数。Occ[c,r]表示在BWT（T）中第r行之前出现字符c的个数。以上面的例子说，C[a]=0,C[c]=3,C[g]=5。Occ[a,5] =2 。记住，是从第0行开始算。明白了C[c],和Occ[c,r]的定义，我们就看一下怎么通过BWT（T）还原出T。

这是LF-Mapping的算法。

这是如何从BWT（T）还原为T的UNPERMUTE algorithm。

这还是上面那个例子举例，描述算法的过程。自己理解一下这个算法，还是不难懂的。

BWT（T）用处之二：可以进行read的精确匹配，EXACTMATCH。

比如，我们在“acaacg”中查找“aac”，就可以用到下面的算法EXACTMATCH（）

LFC（r,c）是对LF（r）的改进，在LF（r）中，第r行对应的只能是字符c，而LFC（r,c）中可以是任意的c。稍微有些差别，不是很重要。按前面的理解也能理解透。

查找的过程就是从后向前找，不断缩小[sp,ep]这个可能匹配的区间，直到找到read的第一个字符a结束。这个[sp,ep]区间就是精确匹配的区间。

2、FM-index的第二部分，checkpoint是个什么东西？

通过上面的介绍，我们知道了什么是BWT（T），什么是LF-Mapping，怎么用BWT（T）还原一个T，怎么查找一个read。

在上面的算法中，无论是EXACTMATCH（）还是UNPERMUTE algorithm，我们都反复的调用OCC（r,c），如果每次调用都从前往后查一遍的话，无疑太浪费时间，幸好，这个OCC（r,c）是固定不变的，所以我们可以把OCC（r,c）提前写好，存在内存中以供查找。但是，对于人类基因组来说，3G的长度，每一行都要有OCC（r,A、C、G、T）无疑开销是十分大的。一个OCC（r,A）用4B，一行就用16B，3G行。。。。可想而知。但是，我们根本没有必要每行都设置四个OCC，只需要隔上一段距离存一行OCC（r,A、C、G、T）即可。这就是一个checkpoint，如果我们查的时候没碰到checkpoint，只需从上一个checkpoint向下捋一遍即可。实际上，在bowtie里每个448行设置一个checkpoint。

3、FM-index第三部分，怎么还有一个SA[]是干什么的？

大家想想，我们讲到这里，知道了怎么查找一个read，但是mapping的终极目标是把read映射到参考基因组上，告诉你这个read在参考基因组上的什么位置。光靠上面这些还是不够的。但是也差不远了。我们只需在下面这个数组第一列前添加一个数组，SA[r]，记录第r行在参考基因组中是什么位置 即可。这个过程实际上是在BWT中实现的。

上面这个例子实际上是BWT的过程，虚线的那个数组就是SA[r]。

对于3G长的串，我们要存SA[]，就得需要12G，这无疑也是负担不起的。能不能利用之前checkpoint的方法，只存放一部分呢？答案是可以的。实际上，在bowtie中每隔36行存一个SA[],因此我叫它简化了的SA[]。

还是之前那个例子，我们找到了aac在哪一行，想通过SA[r]映射到原串上，但恰好没存这里的SA[r]怎么办？那我们就继续进行LF（r）操作，知道碰到了存放SA的哪一行，比如我们又走了两步，才找到带SA的一行i，SA[i]=0,那我们需要的SA[r]=SA[i]+2=2。

OK，写到这里算是完全介绍完了bowtie里的FM-index是怎么个东东。就是一个BWT[T],一个checkpoint data，一个简化了的SA。我们看看要对人类基因组建立FM-index需要多大空间呢。一个碱基可以用2bit表示，那么3G长的BWT大概需要680M空间。Checkpoint data需要BWT[T]的14%左右，简化SA需要BWT[T]的50%左右。加起来的话，建立一个FM-index需要1.1G，实际上在bowtie里存放了两条FM-index，因此需要2.2G左右的空间。我们现在的机器大概都是4G左右，所以完全可以在我们的电脑里运行bowtie。

（author QQ：大海浪涛）