谈谈对BWT算法的理解
最近对分析流程里用的bwa和mutect2的原理感兴趣,打算稍微了解下。说起来本科时候也学过一些经典算法,印象最深刻的就是NW全局比对和SW局部比对了。当时老师是正好抽到我的代码作业,好像我是有个回溯没整出来,最后得个70分。就是工作后发现这东西压根用不上。(而且讲真,如果精于算法和统计的话,生统和算法工程师,待遇和发展上都比生信强。)总之业余选手一名,如有谬误还请指出。
BWT算法的作用,是使字符串中相同字母临近,便于后续压缩。
而平时用到最多的bwa-mem比对,其实是在逆向解码的过程。
首先,假定要处理序列为 ACAACG,将序列末尾补上$
即:ACAACG$
迭代将字符串的首字符,移到末尾,直到$到达首位。最终得到原始的矩阵如下(也有说将末尾字符移到首位的,这个其实不影响)。
image.png
现在,我们先不管原始矩阵排序的事,取出原始矩阵中的首列(F列)和尾列(L列),可以得出以下信息:
第一行:$后的字符为A,这个A是所在F列中的第一个A
第二行: A在字符C前,A是所在L列中的第一个A,C是所在F列中的第一个C
第三行:C在字符A前,C是所在L列中的第一个C,A是所在F列中的第二个A,
第四行:A在字符A前,一个A是所在L列中的第二个A,另一个A是所在F列中的第三个A
第五行:A在字符C前,A是所在L列中的第三个A,C是所在F列中第二个C
第六行:C在字符G前,C是所在L列中第二个C,G是所在F列中第一个G
第七行:G在字符$前,G是所在L列中第一个G
而显然F列与L列在删去$后完全相同,也就是说L列中的第n个A/T/C/G就是F列中的第n个A/T/C/G。
显然,根据这些信息,可以反推出原始序列为 ACAACG。
现在,我们将原始矩阵每行按字典排序
image.png
同理,当然也可以推出原始序列。
已知对每行来说,L列的字符都在F列的前面。则从F列的结尾字符$开始回溯(见下图绿框)。
得到 $ ->G->C->A->A->C->A->$,即ACAACG
只是这时候有一个我也不了解的问题,显然上述成立要求是——在字典排序后,L列中的第n个A/T/C/G就是F列中的第n个A/T/C/G,但这一论点,要如何证明成立?只能说,这可能是字典排序的特性?不确定的事我就不乱说了哎。