BWT算法解码

博文复制了原博客，做了少许修改

将原始序列经过BWT转换后， 可以更方便的进行压缩；而且BWT转换是一个可逆的转换，能够根据转换后的序列还原出原始序列；

BWT转换首先将序列进行在序列的末尾插入一个字符，并且规定按照字典序的排序的话， 这个字符小于序列中的任意字符：

比如原始序列:

acaacg

首先在末尾添加一个$符号，变成 acaacg$;

对这个序列进行循环右移，形成一组序列，；

1) acaacg$

2) $acaacg

3) g$acaac

4) cg$acaa

5) acg$aca

6) aacg$ac

7) caacg$a

然后对这7个序列按照字典序进行排序：

　First     Last

1) $acaacg

2) aacg$ac

3) acaacg$

4) acg$aca

5) caacg$a

6) cg$acaa

7) g$acaac

对于上面的这个序列矩阵，有以下几个性质：

1）last列的第一个字母就是原始序列的最后一个字符；(因为$小于任意字符)

2）first 列的每一个字母是对应的Last列的后一个字符； (因为循环右移)

3）字符x在last列中第i次出现对应first 列中的第i次出现的那个字符，对应的L列便是该字符的上一个字符；比如last列中第5行的a是第2次出现的a,它对应的就是First列中的第3行的a(此处的a在First列中也是第2次出现！），对应到L列中$，因此该a是第一个字符；

利用上面的这3个性质就可以根据First 列和Last 列中的内容，推测出原始序列， 下面是bowtie论文中给出的示意图：

L列第一个字符g为原字符串最后一个字符，依次向前查询，在F列中找到对应字符，就可以获得上一个字符；

利用这个性质进行精确匹配：

比如查询字符串为aac, 参考字符串为acaacg;

1）首先从查询字符串的最后一个字符 c  开始，c 在First 列中出现了两次；这两个c 对应的Last 列中的字符都为a ， 与查询序列相同；

2）检测这两个a 对应的前一个字符， 一个对应$; 一个对应a, 很明显对应a的就是我们想要找到的， 这样查询序列与参考序列就精确匹配上了；

非精确匹配(允许错配)：

在实际的分析中， 由于测序错误或者snp等；要求匹配是必须允许错配和gap;

1) 测序错误

参考基因组序列为ATCCG, 而我们测出来的为AACCG, 由于测序错误， 将第二个碱基T测成了A, 在mapping的过程中，程序应该可以校正这种错误， 事先规定错配的最大碱基数， 这样可以有效的校正一些个别的测序错误；

2) snp

基因组中经常会发生snp,即单核苷酸多态性，对应的碱基发生了突变，可能是A变成了T(转换), 或者A变成了C(颠换),这样的位点在比对的过程中，如果不允许错配，就导致这样的序列比对不到参考基因组上， 这也是不合理的；广义上的snp还包括插入和缺失，这样在比对的过程中，必须允许出现gap,这样才能将序列正确的比对到参考基因组上， 在下游的call snp的分析过程中才可以正确的识别snp位点， 只有bowtie2 支持gappe d 比对， bowtie1 不支持；

bowtie也能够进行非精确的匹配：

在实际比对过程中， 要能够给定比对上的序列在参考基因组上的位置， 这样需要我们存储位置信息，

以上面的BWT转换形成的矩阵为例，需要建立一个数组， 这个数组保存的是对应的字符在First 中第一次出现的位置和该字符在Last列中是第几次出现；

a[0] = (5,1); (g 在 First 列中第一次出现在第5行， 在Last列中是第一次出现的g, 5+ 1得到其在参考序列上的位置为6，)

a[1] = (4,1); 

a[2] = (0,1);

a[3] = (1,1);

a[4] = (1,2);

a[5] = ();

a[6] = ();