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
S_i
Si 代表第
i
i
i 个核苷酸,
1
≤
i
≤
N
1≤i≤N
1≤i≤N ,
N
N
N 为核苷酸序列的长度。
S
i
j
S_{ij}
Sij 代表从
S
i
S_i
Si 到
S
j
S_j
Sj 的子序列 ,
1
≤
i
<
j
≤
N
1≤i<j≤N
1≤i<j≤N。图中画出了核苷酸之间的内部结构,只有G-C、A-U、G-U才能相连。需注意两核苷酸之间的键不能相交。
Background
结构的自由能与键之间的区域,即特征有关。定义 E ( k ) E(k) E(k) 为该特征 k k k 所对应的自由能。结构的自由能为所有特征自由能之和,作者的目的为找到自由能最少的结构。将不可能出现的特征所对应的能量参数设为无穷大,便可避免结构中出现这些特征。
Method
定义
W
(
i
,
j
)
W(i, j)
W(i,j) 为子序列
S
i
j
S_{ij}
Sij 所有可能结构的最小自由能。当
S
i
S_i
Si 和
S
j
S_j
Sj 为碱基对时,定义
V
(
i
,
j
)
V(i, j)
V(i,j) 表示
S
i
j
S_{ij}
Sij 所有可能结构的最小自由能。当
S
i
S_i
Si 和
S
j
S_j
Sj 不是碱基对时
V
(
i
,
j
)
=
∞
V(i, j) = \infty
V(i,j)=∞ 。作者递归的计算
V
(
i
,
j
)
V(i, j)
V(i,j) 和
W
(
i
,
j
)
W(i, j)
W(i,j) ,先从五核苷酸序列开始,直到越来越长的核苷酸序列。因为五核苷酸序列没有稳定的结构,先验知识, 所以当
j
−
1
=
4
j-1=4
j−1=4 时
W
(
i
,
j
)
=
0
W(i, j) = 0
W(i,j)=0 。
定义
F
H
(
i
,
j
)
FH(i, j)
FH(i,j) 表示由
S
i
S_i
Si 、
S
j
S_j
Sj 形成的键组成的hairpin环特征,定义
F
L
(
i
,
j
,
i
′
,
j
′
)
FL(i, j, i', j')
FL(i,j,i′,j′) 表示由
S
i
S_i
Si 、
S
j
S_j
Sj 形成的键和
S
i
′
S_i'
Si′ 、
S
j
′
S_j'
Sj′ 形成的键组成的特征。
当特征仅含有一条键时,如图2-A1所示,
V
(
i
,
j
)
=
E
(
F
H
(
i
,
j
)
)
V(i, j) = E(FH(i, j))
V(i,j)=E(FH(i,j)) 。
当特征包含两条键时,如图2-A2所示,
V
(
i
,
j
)
=
E
(
F
L
(
i
,
j
,
i
′
,
j
′
)
)
+
V
(
i
′
,
j
′
)
V(i, j) = E(FL(i, j, i', j')) + V(i', j')
V(i,j)=E(FL(i,j,i′,j′))+V(i′,j′) ,其中
i
<
i
′
<
j
′
<
j
i<i'<j'<j
i<i′<j′<j 。
当特征包含两条以上的键时,如图2-A2所示(bifurcation环),
V
(
i
,
j
)
=
W
(
i
+
1
,
i
′
)
+
W
(
i
′
+
1
,
j
−
1
)
V(i, j) = W(i+1, i') + W(i'+1, j-1)
V(i,j)=W(i+1,i′)+W(i′+1,j−1) ,其中
i
+
1
<
i
′
<
j
−
2
i+1<i'<j-2
i+1<i′<j−2 。bifurcation环的计算使用interior环的计算方法近似
V
(
i
,
j
)
V(i, j)
V(i,j) 表示
S
i
j
S_{ij}
Sij 所有可能结构的最小自由能,表达式如下
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_i
Si 和
S
j
S_j
Sj 至少有一个不在键上时,如图2-B1所示,
W
(
i
,
j
)
=
W
(
i
+
1
,
j
)
W(i, j) = W(i+1, j)
W(i,j)=W(i+1,j) 或
W
(
i
,
j
)
=
W
(
i
,
j
−
1
)
W(i, j) = W(i, j-1)
W(i,j)=W(i,j−1) 。
当
S
i
S_i
Si 和
S
j
S_j
Sj 在同一条键上时,如图2-B2所示,
W
(
i
,
j
)
=
V
(
i
,
j
)
W(i, j) = V(i, j)
W(i,j)=V(i,j) 。
当
S
i
S_i
Si 和
S
j
S_j
Sj 在不同键的两侧时,如图2-B3所示,
W
(
i
,
j
)
=
W
(
i
,
i
′
)
+
W
(
i
′
+
1
,
j
)
=
W
(
i
,
j
′
−
1
)
+
W
(
j
′
,
j
)
W(i, j) = W(i, i') + W(i'+1, j) = W(i, j'-1) + W(j', j)
W(i,j)=W(i,i′)+W(i′+1,j)=W(i,j′−1)+W(j′,j) ,其中
i
<
i
′
<
j
′
<
j
i<i'<j'<j
i<i′<j′<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
有太多不同结构的特征具有相似的自由能,作者认为以后发现的一些规则和信息能够融合进这个模型中,帮助模型变得更为准确。