【基因序列比对】Needleman-Wunsch算法

前言

最近加入了一个研究基因大数据的项目组，要学习基因相关知识QAQ，所以开始一个新系列梳理知识和激励自己学习。

本系列参考了https://zhuanlan.zhihu.com/p/54142276

Needleman-Wunsch算法

Needleman-Wunsch算法是基于动态规划的算法，记不清楚的可以看下我之前的动态规划的文章，简单总结就是将一个问题分解为很多相互联系并且逐渐扩大的小问题，依次解决使问题规模逐渐扩大，最终获得问题的解。

DNA有三种情况可能导致两个序列不同：

1. SNP单核酸多态性，就是碱基被替换了，出现频率最高（CGTT -> CGTA）
2. INSERT，就是多复制了一个（CGTT -> CGTTT）
3. DELETION，就是少复制了一个（CGTT -> CGT）

那么两个序列的碱基一一比对时共有三种可能的情况

1. MATCH：上下匹配
2. MISMATCH：出现SNP，上下不匹配
3. GAP：出现INSERT或DELETION导致一个序列为空缺。

我们在比对时需要一定的准则来评判不同错误产生的比对损失，因为在测出两个序列的时候我们不知道引起发错误的原因。

假定MATCH得1分，MISMATCH或GAP得-1分。我们需要知道的就是两个序列在哪种排列对应时能获得最高分。我们当然可以使用穷举法获得最优解，但是序列长的话费时费力，那么分析情况会发现，使用动态规划妙不可言！

举个栗子

假定要比对的两个序列：
ATGCATG
AACCGTC

步骤1：初始化得分表
由于Needleman-Wunsch算法是基于动态规划算法思想的，所以他必有一个得分表可用于记录和回溯。
第一行和第一列为初始的得分，因为相邻构成一个GAP，例如第一行第三列的-2：-AT与—有两个GAP，故为-2。

步骤2：填表
用从左上到右下的顺序计算每个位点的得分，每个位点的得分与左上、左、上方向的得分有关，计算规则为移动前位点分数和移动造成的得分相加，选取三个方向得分的最高分（即最优解）为当前位点得分。依此类推填完所有表格。

以第一步为例：
左上：
原分数：0分
移动造成分数：

• -> -A
• -> -A
  为MATCH：得1分
  总分数：0 + 1 = 1分

上：
原分数：-1分
移动造成分数：
-A -> -A-
– -> --A
造成GAP：-1分

总分数：-1 + (-1) = -2分

左：
原分数：-1分
移动造成分数：
– -> --A
-A -> -A-
造成GAP：-1分

总分数：-1 + (-1) = -2分

选择其中最高分：左上方向的1分，填入表中即可

依此类推填满所有表格：

回溯找回原来的路径：
从右下方开始，如果选值不为左上方，则在相应方向引入一个GAP，使最终路径为一个对角线即可得到最终的最优对应序列。
按照上法可得到三个比对得分都为0的结果（最优解）：
原序列
ATGCATG
AACCGTC
最优序列
ATGCA-TG
A-ACCGTC

ATG-CATG
A-ACCGTC

ATGC-ATG
A-ACCGTC