今天看了下bowtie 的论文, 里面描述了BWT转换的过程和bowtie的比对算法;
NGS测序数据的数据量非常大, 为了更快的处理, 通常需要对数据进行压缩;而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个字符;比如last列中第5行的a是第2次出现的a,它对应的就是First列中的第3行的a;
利用上面的这3个性质就可以根据First 列和Last 列中的内容,推测出原始序列, 下面是bowtie论文中给出的示意图:
利用这个性质进行精确匹配:
比如查询字符串为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] = ();