全局比对与动态规划

前言

学过生信的肯定知道Needleman–Wunsch算法和Smith–Waterman 算法,一个用来进行全局比对,一个用来进行局部比对。单纯看算法抽象后的算法公式的话也不复杂,用两个短序列来套公式计算的话,画过箭头的同学都知道不难。本文简单讲述Needleman–Wunsch算法是如何利用动态规划来寻找两序列最大相似度。

Python代码实现:Needleman–Wunsch 算法代码

动态规划

动态规划是数学上的优化方法,也是计算机编程的一种方法。对于一名生物狗,从数学层面去看动态规划,就太过硬核了。我们就从计算机编程角度理解动态规划就好了。
动态规划是一种解决问题的方法,将一个复杂大问题拆分为几个小的问题,从几个小问题的最优解推出大问题的最优解。当然这不是说将一个大问题拆分了就能用动态规划了。同时问题也需要满足两个性质,这两个性质后面再讲。先回到我的主要探讨问题全局比对上:对于两条蛋白质或者核酸序列来说,寻找它们的最大相似度,如何拆分成几个小问题呢?(提示:矩阵、算法公式)

拆分问题

就这样凭空想,如何拆分,我觉得我是想破脑袋也想不了T_T. 所以我们按照,学全局比对的套路手动来一遍序列比对操作,从中看看如何拆分问题。

为方便计算,得分矩阵就弄个简单点的:

N A T G C -
A 10 5 5 5 -5
T 5 10 5 5 -5
G 5 5 10 5 -5
C 5 5 5 10 -5
- -5 -5 -5 -5 0

也就说完全匹配得10分,不全得5分,匹配到空位-5分,空位对空位一开始用来初始化时,才有意义。在之后比对中空位对空位没有意义。
先给出公式(d为gap罚分,即碱基比对空格的得分):


来自高歌老师PPT

现假设我们有两条序列Seq1: CATCCCCCAAA 和 Seq2: CAAAA。长度分别为11,5。根据它们的长度画一个合适的二维矩阵,填充好初始数据:


image.png

加粗的分别是 Seq1,Seq2,其中 -代表空格。一开始空格对空格,我们初始化,记F(0,0)=0,同时记空格分别对应的行和列是第0行,第0列。然后第1行,第1列分别是两条首字母对应的行列。用(i,j)代表(行,列),i,j能有多大,取决于序列对的长度。s(x_i,y_i)的分数是一条序列中第i个碱基和另一条序列中第j个碱基比对的分数。放在矩阵中就是了第i行的碱基和第j列的碱基比对结果。

按照F(i,j)函数公式,F(i,j)的值由左上的值(F(i-1,j-1)),上方值(F(i-1,j))和左方值(F(i,j-1))确定。
(注:这里的左上,上,左是一个相对位置,在当前设定矩阵坐标系中才有意义,这里是将(0,0)位置放在了左上方,向右和向下延伸。但其实,你也可以将出放在其它位置,比如右下方)

极端比对

我们看下第0行,它们只能由左方的值得到。我们算几个值:
F(0,1)=F(0,1-1)-5=0-5
F(0,2)=F(0,2-1)-5=-5-5=-10
F(0,3)=F(0,3-1)-5=-10-5=-15
第0列,只能由上方的值得到,计算同理,不再演示了。
因为第0行和0列,很容易就得出来了,所以一般画矩阵时,就顺便一起画出来了。

你应该知道,往右或往下延申,则比对到了空格。那这是为啥?
就拿第0行的比对情况举个例子,第0行中,i=0已经固定了,它代表一条序列(后面称为S1序列)的第0元素。j=0,1,2,3...,分别代表另一条序列(后面称为S2序列)第0,1,2,3...元素。序列第一个元素是序列的第一个碱基。第0个是啥玩意,第0个啥也不是,那我们就把它当做空格吧。

  1. 序列比对从两条序列的第(0,0)元素开始,两个空格比对,得分为0。
  2. i=0固定,一直表示S1序列第0个元素,但是j是一直增加的,所以j=1时,就是用S2序列的的第j个元素和S1序列的第i=0个元素比对,但是S1序列第i=0元素已经和S2序列的第j=0元素比对过了,你不能重复比对啊,你也不能私自和S1序列的i=1元素比对,所以没办法,就只能和一个空格比对,产生空位罚分,在之前得分基础上上加上空位罚分。
  3. S2序列的第2,3,4...元素都是如此,只能比对空格,产生空位罚分。

比对往右延申的过程就是行数i固定了,没有增加,但是列数j增加,又两条序列的中的元素不能重复比对,也不能比对序列其它的元素,所以就只能比对空格了。(向下延申同理。)
我们看两个极端的例子:

image.png

##黑色路线
CATCCCCCAAA-----
-----------CAAAA
##红色路线
-----CATCCCCCAAA
CAAAA-----------

你们说这第0行,0列的有啥用处,放这么多负分而产生极端比对。要排除这些不合适比对结果也要花费许多资源。

比对开始

上面进行了一次小探讨。比对直接向右或向下会产生gap。因为行列有一方被固定,没法提供可比对元素。如果要正确比对的话,要使得ij同时增加,即向右下方延申。除了第0行和0列的元素,只有单一来源,矩阵其它的元素都有三个来源,即上方,左方,左上。如何确定F(i,j)的值,按照公式,选择其中最大的值。F(m,n)(m,n为两条序列的长度)为最终序列比对的相似得分。
既然我们想使F(m,n)最大,肯定不能过多的产生空位罚分吧。那么我们从一开始就一直向右下方延申,直至不能继续,能到得到最大的相似得分呢?我们试一下吧:

image.png

我们用箭头表示从哪个方向得到的最大值,两个箭头则可以从两个方向得到。从图中红色箭头的路线可以看出来,一直选择往右下延申的比对,似乎可以达到最大值。而且再比对三对元素就可以到达终点了。现在看红色箭头,似乎是全场全场最靓的仔。然而最终结果如何?我们继续把剩余几个各自补全看看。
image.png

现在再看全场最靓的仔似乎被绿了??而且是一个绿一个的。简直就是螳螂捕蝉,黄雀在后,然后渔翁得利。。。

贪心算法:红色路线选择的策略,从一开始,三个方向的选择,每一步,它都选择三个方向中得分最大的那个,但最后结果却不是最好的。这种在每一步都做出当前看来是最好选择的方法,而不考虑整体最优的方法,叫做贪心算法。这种方法在这里得不到最优结果。
穷举搜索:这是一种很直观的方法,就是把所有可能的情况都列举出来,从中找最优的。但是用于序列比对是是非常差的一种方法。从矩阵左上方到达矩阵右下方,有多少种路线?矩阵每一个小格有三种可能路线,矩阵大小m*n,即有3^{mn}种路线。其中m,n稍微取大一点的值,数值就爆炸了。完全不可取。

emm,写到这里,还没讲到全局比对是如何拆分问题的,现在我们就开始讨论这个问题吧。。
首先,在矩阵中当比对到第(m,n)格时,代表着比对终止。因为两条序列的每个元素都参加了比对过程,无论序列中的元素是比对到了空格,还是另一条序列的某个元素。因为序列中的元素都比对完了,在比对下去就只有空格对空格了,没什么意义。所以(m,n)是终点。
我们就从终点开始探讨如何分割大问题为小问题吧。
F(m,n)=max\begin{cases}F(m,n)+s(x_i,y_j)\\F(m-1,n)+d\\F(m,n-1)+d\end{cases}
其中s(x_i,y_j)d是知道的,那问题F(m,n)则可以通过求出F(m-1,n-1),F(m,n-1)F(m-1,n)的分数得到。即一个大问题,分成了3个子问题了。三个子问题也不知道分数,就继续分子问题呗。如图:

image.png

分治算法:分治的思想也是将一个问题分为几个问题来解决。这与动态规划的分离子问题有啥区别?从上面的图我们可以看出来由F(m,n)分割的子问题们是存在重复问题的。分治分离的子问题一般独立的不会相互重叠,而动态规划的子问题则一般发生重叠现象。

两个性质

之前空着没讲的两个性质,现在来讲一下。一个是最优子结构,一个是无后效性
最优子结构:如果一个问题的最优解,可以将其分为几个子问题,然后能从子问题的最优解推出大问题的最优解。就称其有最优子结构性质。这一性质很容易从函数公式中看出来。
无后效性:给定某一阶段的状态,这一阶段以后的发展过程不受这阶段以前各段状态的影响(来自阮行止的回答)这话是什么意思呢,以我们序列比对为例。比如:当你确定了前面某个阶段的比对,比对情况如下面的比对方案的前面一部分,这时比对分数已经确定了。无后效性就是前面的比对情况对后面的比对结果的分数不会有任何影响。前面部分的比对方案的任意变动也不会影响后面部分的结果。

##1
CAT  CCCCC-----AAA
CA-  ----------AAA
##2
CAT  CCCCCAAA-----
CA-  AA----------A
##3
CAT  CCCCCAAA-----
CA-  ----------AAA

无后效性,可以让我们放心大胆的在矩阵每个格子存放最优的结果,舍弃许多不是最优的分数。因为不存在前面阶段的最差比对方案,还能与后面阶段的最差方案组合成为全局最优比对方案。

复杂度

我们来看一下Needleman–Wunsch算法的复杂度。从二维矩阵中,可以看到我们需要计算每一格的所占最优分值,同时也为了之后的计算方便也得储存每一个的分值。故时间复杂度和空间复杂度都是O(mn)

其它

为了理解动态规划,我在知乎上这个问题什么是动态规划(Dynamic Programming)?动态规划的意义是什么?的答案下看了许多回答。许多答主都给了自己的理解。他们的回答有相似之处,也有不同。不能说谁对谁错(小声bb:也不敢说谁对谁错)
就全局比对我也谈谈我的看法:
动态规划和其它方法有许多相似的地方。要是将其它方法的定义定广一点,也能把动态规划说成其它方法。比如:
穷举搜索: 之前上文曾提到过,不过那时我们说要穷举的对象是每一个比对方案。现在我们换一个观点想一想:穷举每一个矩阵中的格子,即穷举每一个最优的F(i,j)。这里我们的对象就从比对路径而换成了F(i,j)(序列比对到第(i,j)的最大相似度)
贪心算法: 之前提到的贪心算法,我们狭隘的将最优,认为是向下、向右或向右下而能得到的分数。现在我们将最优看成相邻F(i,j)中最大的那个。即按照这个贪心选取,它会可能不断的对比相邻的F(i,j)而选择最优。直至,到最后剩下三条路径到达F(m,n),选中最优的那个。
总之,不管黑猫白猫,能帮助解决问题就是好猫!选择你认为最好的一种理解。在以后多次遇到动态规划后,在慢慢完善你的认知。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 216,692评论 6 501
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,482评论 3 392
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 162,995评论 0 353
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,223评论 1 292
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,245评论 6 388
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,208评论 1 299
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,091评论 3 418
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,929评论 0 274
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,346评论 1 311
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,570评论 2 333
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,739评论 1 348
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,437评论 5 344
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,037评论 3 326
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,677评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,833评论 1 269
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,760评论 2 369
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,647评论 2 354

推荐阅读更多精彩内容