Dynamic programming algorithm

Optimal computer folding of large RNA sequences using thermodynamics and auxiliary information

Year: 1981
Authors: Michael Zuker and Patrick Stiegler
Journal Name: Nucleic Acids Research

Motivation

将动态规划算法与热动力学数据结合。

Research Objective

基于动态规划算法提出一种更有效、快速、并且能折叠更长核苷酸序列的模型。

Dataset

定义 $S_i$ 代表第 $i$ 个核苷酸， $1\le i\le N$ ，$N$ 为核苷酸序列的长度。 $S_{ij}$ 代表从 $S_i$ 到 $S_j$ 的子序列 ，$1\le i<j\le N$。图中画出了核苷酸之间的内部结构，只有G-C、A-U、G-U才能相连。需注意两核苷酸之间的键不能相交。

Background

结构的自由能与键之间的区域，即特征有关。定义 $E(k)$ 为该特征 $k$ 所对应的自由能。结构的自由能为所有特征自由能之和，作者的目的为找到自由能最少的结构。将不可能出现的特征所对应的能量参数设为无穷大，便可避免结构中出现这些特征。

Method

定义 $W(i,j)$ 为子序列 $S_{ij}$ 所有可能结构的最小自由能。当 $S_i$ 和 $S_j$ 为碱基对时，定义 $V(i,j)$ 表示 $S_{ij}$ 所有可能结构的最小自由能。当 $S_i$ 和 $S_j$ 不是碱基对时 $V(i,j)=\infty$ 。作者递归的计算 $V(i,j)$ 和 $W(i,j)$ ，先从五核苷酸序列开始，直到越来越长的核苷酸序列。因为五核苷酸序列没有稳定的结构，先验知识， 所以当 $j-1=4$ 时 $W(i,j)=0$ 。
定义 $FH(i,j)$ 表示由 $S_i$ 、 $S_j$ 形成的键组成的hairpin环特征，定义 $FL(i,j,i^{\prime },j^{\prime })$ 表示由 $S_i$ 、 $S_j$ 形成的键和 $S_i^{\prime }$ 、 $S_j^{\prime }$ 形成的键组成的特征。
当特征仅含有一条键时，如图2-A1所示， $V(i,j)=E(FH(i,j))$ 。
当特征包含两条键时，如图2-A2所示， $V(i,j)=E(FL(i,j,i^{\prime },j^{\prime }))+V(i^{\prime },j^{\prime })$ ，其中 $i<i^{\prime }<j^{\prime }<j$ 。
当特征包含两条以上的键时，如图2-A2所示（bifurcation环）， $V(i,j)=W(i+1,i^{\prime })+W(i^{\prime }+1,j-1)$ ，其中 $i+1<i^{\prime }<j-2$ 。bifurcation环的计算使用interior环的计算方法近似
$V(i,j)$ 表示 $S_{ij}$ 所有可能结构的最小自由能，表达式如下
V ( i , j ) = m i n { E 1 , E 2 , E 2 } E 1 = E ( F H ( i , j ) ) E 2 = min ⁡ i ′ , j ′ { E ( F L ( i , j , i ′ , j ′ ) ) + V ( i ′ , j ′ ) } E 3 = min ⁡ i ′ { W ( i + 1 , i ′ ) + W ( i ′ + 1 , j − 1 ) } V(i, j) = min\{E1, E2, E2\} \\ E1 = E(FH(i, j)) \\ E2 = \min_{i', j'}\{E(FL(i, j, i', j')) + V(i', j')\} \\ E3 = \min_{i'}\{W(i+1, i') + W(i'+1, j-1)\} V(i,j)=min{E1,E2,E2}E1=E(FH(i,j))E2=i′,j′min​{E(FL(i,j,i′,j′))+V(i′,j′)}E3=i′min​{W(i+1,i′)+W(i′+1,j−1)}
当 $S_i$ 和 $S_j$ 至少有一个不在键上时，如图2-B1所示， $W(i,j)=W(i+1,j)$ 或 $W(i,j)=W(i,j-1)$ 。
当 $S_i$ 和 $S_j$ 在同一条键上时，如图2-B2所示， $W(i,j)=V(i,j)$ 。
当 $S_i$ 和 $S_j$ 在不同键的两侧时，如图2-B3所示， $W(i,j)=W(i,i^{\prime })+W(i^{\prime }+1,j)=W(i,j^{\prime }-1)+W(j^{\prime },j)$ ，其中 $i<i^{\prime }<j^{\prime }<j$ 。
W ( i , j ) = m i n { W ( i + 1 , j ) , W ( i , j − 1 ) , V ( i , j ) , E 4 } E 4 = min ⁡ i ′ { W ( i , i ′ ) + W ( i ′ + 1 , j ) } W(i, j) = min\{W(i+1, j), W(i, j-1), V(i, j), E4\} \\ E4 = \min_{i'}\{ W(i, i') + W(i'+1, j) \} W(i,j)=min{W(i+1,j),W(i,j−1),V(i,j),E4}E4=i′min​{W(i,i′)+W(i′+1,j)}
递归算法将核苷酸一个一个的加到序列中，计算每轮的最优结构。本人认为可能陷入局部最优。

Future Work

有太多不同结构的特征具有相似的自由能，作者认为以后发现的一些规则和信息能够融合进这个模型中，帮助模型变得更为准确。