最近对分析流程里用的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
image.png
只是这时候有一个我也不了解的问题，显然上述成立要求是——在字典排序后，L列中的第n个A/T/C/G就是F列中的第n个A/T/C/G，但这一论点，要如何证明成立？只能说，这可能是字典排序的特性？不确定的事我就不乱说了哎。
网友评论