在50亿年后太阳变成红巨星之前能精确比对多少分子序列?

by 林溪|D-Horse

最近在构思一篇关于“计算人类末日”的文章,搜集了一些相关的资料,突然联想到了标题中的这道曾经的作业题,遂随手写下此文。。

(作者不是标题党,不过标题中的这个问题我还是要放到本文的最后再来回答。)

当分子生物学中唯一的法则“中心法则”饱受系统生物学家的攻击和批评[1]的同时,由基因组学发展出来的另一个“中心法则”——序列决定结构,结构决定功能——则显得更加无懈可击。而生物信息学也正是基于这一新的法则来分析海量的生物学数据的。由此,序列分析也成为了这一新学科的基础。这里的“序列”主要指的是生物大分子中的DNA、RNA和蛋白质。

序列分析中的主要工作之一是序列比对。通过检验不同序列间的相似性程度,可以对它们之间的同源性和进化关系进行判断。我们现在看到的那些进化树,虽然在一些细枝末节上还存在很大争议,但它们大多都是通过这种分子水平上的序列比对而构建起来的,比以往的博物学家们只能在形态学和解剖学上进行比较要精确得多了。

通过测序,我们可以把这些线性的生物大分子变成由特定字符集组成的字符串,存储到熟悉它们的计算机中。如DNA的字符集是{A,T,C,G},RNA是{A,U,C,G},蛋白质则是那20种基本氨基酸的单字母缩写。在计算机中用处理字符串的方法来分析这些序列,比用纯生物学的方法要方便、快捷、准确多了。

那么到底要如何比对这些序列呢?这里就不得不要提到一个算法——动态规划(Dynamic Programming).讲解“动态规划”的中文维基页面上用斐波那契数列的例子来解释这个算法是如何工作的。如果没有理解,也没有关系,因为心理学告诉我们,只要能把抽象的问题具体化,其实也没有什么是难理解的。所以这里让我们来考虑一下下面这个更加现实的问题。

假设你要去一个陌生的大城市旅游,但坏消息是你只有不到一天的时间,肯定不能把城市里所有的景点游览完。同时还有个好消息:你有这个城市的地图,可以提前规划路线。作为限定条件,在这个城市里游览的时候,你只能向东或向南前进,不能走回头路。那么你要如何选择路线才能尽可能多地浏览景点呢?这个问题也就是经典的“曼哈顿游客问题”(Manhattan tourist problem).

“曼哈顿游客问题”中的地图。(在这里我们可以忽略那条斜街。)图像来源: An introduction to bioinformatics algorithms

假设我们从地图的西北角进入,终点是地图的东南角。利用动态规划算法求解的主要思路就是计算到达每一个汇点时所能游览景点的最大数,也就是将原问题分割成相同的子问题,再分别对子问题依次求解。将“曼哈顿游客问题”一般化,引用下图来继续说明。图中箭头代表街道和前进的方向,箭头上方的数字代表这条街道上的景点*加权值(比如相对历史博物馆我可能更喜欢自然博物馆,那么就可以设自然博物馆的加权值为5,历史博物馆为1)。圆圈代表两条路的汇点,圈中数字代表之前游览过的景点*加权值的最大值总和。我们在计算每一个汇点的时候,都要考虑能使该汇点达到最大值的路线。将地图计算完成的情况如下方左图所示,粗箭头即代表选择的行进方向,最后通过回溯的方法,一条最优路线就浮现出来了(下方右图所示)。

“曼哈顿游客问题”动态规划算法求解一般过程。图像来源: An introduction to bioinformatics algorithms

 在序列比对中,我们也沿用这一思路,不过要先定义一些与序列分析相关的操作:

与“曼哈顿游客问题”中景点的“加权值”类似,我们还需要一个罚分矩阵:

也就是说对每一次“匹配”操作记0分,对其它操作均记1分。事实上这是一个简化了的罚分模型[2]。

现在我们就可以对序列进行比对了。为简化计算,假设有这样两条“DNA超短序列”:u序列TTCACGv序列TGTACAG.与在“曼哈顿游客问题”中我们需要寻找一条最大值的路径不同的是,在这里我们想寻找的是一条满足最小值的路径,即最小罚分。

算法描述如下:

For 1≤i≤|u|,1≤j≤|v|,

计算如下图所示:

在上图动态规划计算中,一个斜箭头表示一次“匹配”或“错配”操作,一个向上的箭头表示一次u序列对v序列的“插入”操作,一个向左的箭头表示一次u序列对v序列的“缺失”操作。最后值为3,也就是此次序列比对在最大化进行匹配情况下的最小分值。

最后沿箭头方向进行回溯,在本例中可得两个最优解:

以上就是利用动态规划算法对两条序列进行比对的一般求解过程。假设两条序列长度均为n,那么空间复杂度很明显就是O(n^2),时间复杂度,也就是总共需要计算的次数。对于每个中部节点,我们每次要计算3次,一共n^2个中部节点,所以时间复杂度为O(3*n^2),可以简化记为O(n^2).

以上我们考虑的序列比对问题都是“双序列比对”(pairwise alignment),即只同时比对两条序列。那么我们现在在双序列比对的情况下,来回答下标题的那个问题吧。这里还需要补充说明的是,一天计86400秒,一年计365.2425天,假设所有序列的长度均为1000.事实上绕了上面那一大圈,解决这个问题并不复杂。不考虑超级计算机也不考虑并行计算,我们只考虑一台极普通的家用电脑的运行能力,每进行一次双序列比对的时间是2秒,那么一共可计算的序列对数就是:

86400*365.2425*5*10^9/2=7.889238*10^16(对)

10的16次方级,感觉很多不是么,但是本文想讨论的显然不是这个。“双序列比对”固然很有用,但在很多时候,我们更需要的是多序列比对,也就是同时比对3、4条,甚至7、8条序列,那么情况又会如何呢?这就需要把动态规划扩到更高维的空间去了。

高维空间的动态规划。图像来源:Skript von Uni Bielefeld

我们暂时先考虑三维空间的情况。如上图所示,每一个中部节点对应的是7条楞/边,时间复杂度瞬间变为O(7n^3),可见动态规划算法的时间复杂度是随着k序列比对问题中k值的增加而呈指数级增长的。推广至一般化情况,假设是k序列比对问题,即同时比对k条序列,那么时间复杂度就是O[(2^k-1)*n^k],现在再来计算多序列比对在标题中时间限定下的k值(其它条件如前文不变):

(2^k-1)*n^k=86400*365.2425*5*10^9,k=5.21

是的,50亿年,5条序列(如果有人不了解“50亿年后太阳变成红巨星”意味着什么的话可以参考这个维基页面中的描述)。别忘了,这是在精确比对计算的前提条件下得出的结论。当然,我们人类不可能真的耗上这几十亿年去计算5条序列出来,而且事实上在很多基因组级别的序列比对问题中,单序列的长度也远远超过1000,所以这也就成为了生物信息学家们面临的一个挑战:要尽可能地降低算法时间复杂度,同时还要使计算结果具有生物学意义(当然空间复杂度也要想办法降低,但是更为重要的还是使时间复杂度降低。事实上,动态规划算法就是一种用空间换时间的算法,可以将其与递归算法相比较)。

[1]关于“中心法则”的讨论,可以参考这篇Dog eat dogma,译言上有一篇翻译,滑落神坛的中心法则。有人曾经问“中心法则”的提出者佛朗西斯·克里克(即DNA双螺旋结构的发现者之一)为什么直接就称之为“法则”(其实英文为Dogma,有“教条”、“教义”的意思),克里克的回答是因为当时他已经提出了一个假说,叫“序列假说”(Sequence hypothesis),他认为“中心法则”中的内容,即蛋白质的合成路径更为重要,为彰显其重要性,便称其为“法则”。

[2]使序列比对中的罚分模型更有生物学上的意义:两个在序列比对中较为常用的罚分模型,一个是PAM,另一个是BLOSUM.