LinearPartition: linear-time approximation of RNA folding partition function and base-pairing probabilities
Year: 2020
Authors: He Zhang, Liang Zhang, David H. Mathews and Liang Huang
Journal Name: Bioinformatics
Motivation
传统分割方法的复杂度与序列长度呈三次方关系。
Research Objective
设计一种线性时间的算法来近似分割函数和碱基对概率
Dataset
定义 x = x 1 , . . . , x n \textbf{x} = x_1, ..., x_n x=x1,...,xn 为长度为 n 的 RNA 序列, Y ( x ) \mathcal{Y}(\textbf{x}) Y(x) 代表 x \bf{x} x 所有可能的二级结构。
Method
分割函数为
Q
(
x
)
=
∑
y
∈
Y
(
x
)
e
−
Δ
G
∘
(
y
)
R
T
Q(\textbf{x}) = \sum_{\textbf{y} ∈ \mathcal{Y}(\textbf{x})} e^{-\frac{ΔG^{\circ}(y)}{RT}}
Q(x)=y∈Y(x)∑e−RTΔG∘(y)
其中,
Δ
G
∘
(
y
)
ΔG^{\circ}(y)
ΔG∘(y) 为结构
y
y
y 的 Gibb 自由能变化,
R
R
R 是气体常量,
T
T
T 是热动力学温度。
Δ
G
∘
(
y
)
ΔG^{\circ}(y)
ΔG∘(y) 使用逆 Nussinov–Jacobson 能量模型计算。在位置
j
j
j 的未配对碱基的自由能变化表示为
δ
(
x
,
j
)
δ(\textbf{x}, j)
δ(x,j) ,碱基对
(
i
,
j
)
(i, j)
(i,j) 的自由能变化表示为
ξ
(
x
,
i
,
j
)
\xi(\textbf{x}, i, j)
ξ(x,i,j) 。
因此,
Δ
G
∘
(
y
)
ΔG^{\circ}(y)
ΔG∘(y) 表示为
Δ
G
∘
(
y
)
=
∑
j
∈
u
n
p
a
i
r
e
d
(
y
)
δ
(
x
,
j
)
+
∑
(
i
,
j
)
∈
p
a
i
r
s
(
y
)
ξ
(
x
,
i
,
j
)
ΔG^{\circ}(y) = \sum_{j ∈ unpaired(y)} δ(\textbf{x}, j) + \sum_{(i, j) ∈ pairs(y)} \xi(\textbf{x}, i, j)
ΔG∘(y)=j∈unpaired(y)∑δ(x,j)+(i,j)∈pairs(y)∑ξ(x,i,j)
其中,
u
n
p
a
i
r
e
d
(
y
)
unpaired(y)
unpaired(y) 表示
y
y
y 中未配对的碱基集合,
p
a
i
r
s
(
y
)
pairs(y)
pairs(y) 表示
y
y
y 中的碱基对集合。
分割函数可表示为
Q
(
x
)
=
∑
y
∈
Y
(
x
)
(
∏
j
∈
u
n
p
a
i
r
e
d
(
y
)
e
−
δ
(
x
,
j
)
R
T
∏
(
i
,
j
)
∈
p
a
i
r
s
(
y
)
e
−
ξ
(
x
,
i
,
j
)
R
T
)
Q(\textbf{x}) = \sum_{\textbf{y} ∈ \mathcal{Y}(\textbf{x})} (\prod_{j ∈ unpaired(y)} e^{-\frac{δ(\textbf{x}, j)}{RT}} \prod_{(i, j) ∈ pairs(y)} e^{-\frac{\xi(\textbf{x}, i, j)}{RT}})
Q(x)=y∈Y(x)∑(j∈unpaired(y)∏e−RTδ(x,j)(i,j)∈pairs(y)∏e−RTξ(x,i,j))
作者定义
[
i
,
j
]
[i, j]
[i,j] 为子序列
x
i
,
.
.
.
,
x
j
x_i, ..., x_j
xi,...,xj (
[
j
,
j
−
1
]
[j, j-1]
[j,j−1] 为空)。
[
i
,
j
]
[i, j]
[i,j] 对应的分割函数
Q
i
,
j
=
∑
y
∈
Y
(
x
i
,
.
.
.
,
x
j
)
e
−
Δ
G
∘
(
y
)
R
T
Q_{i, j} = \sum_{y ∈ \mathcal{Y}(x_i, ..., x_j)} e^{-\frac{ΔG^{\circ} (y)}{RT}}
Qi,j=∑y∈Y(xi,...,xj)e−RTΔG∘(y)
线性切割算法的伪代码如下所示
- function LinearPartition( x \textbf{x} x, b)
- n = len( x \textbf{x} x )
- 令所有的 Q j , j − 1 = 1 Q_{j, j-1} = 1 Qj,j−1=1
- for j in (1, …, n) do
- for each [ i , j − 1 ] [i, j-1] [i,j−1] in Q Q Q do
- Q i , j + = Q i , j − 1 ⋅ e − δ ( x , j ) R T Q_{i, j} += Q_{i, j-1} · e^{-\frac{δ(\textbf{x}, j)}{RT}} Qi,j+=Qi,j−1⋅e−RTδ(x,j)
- if x i − 1 x j x_{i-1}x_j xi−1xj in {AU, UA, CG, GC, GU, UG} then
- for each [ k , i − 2 ] [k, i-2] [k,i−2] in Q Q Q do
- Q k , j + = Q k , i − 2 ⋅ Q i , j − 1 ⋅ e − ξ ( x , i − 1 , j ) R T Q_{k, j} += Q_{k, i-2} · Q_{i, j-1} · e^{-\frac{\xi(\textbf{x}, i-1, j)}{RT}} Qk,j+=Qk,i−2⋅Qi,j−1⋅e−RTξ(x,i−1,j)
- BeamPrune( Q Q Q, j, b)
- return Q Q Q
在 第 5-6 行,通过添加碱基
x
j
x_j
xj 将
[
i
,
j
−
1
]
[i, j-1]
[i,j−1] 扩展为
[
i
,
j
]
[i, j]
[i,j] 。在 8-9 行,如果
x
i
−
1
x_{i-1}
xi−1 和
x
j
x_j
xj 为碱基对,就结合
[
i
,
j
−
1
]
[i, j-1]
[i,j−1] 和左侧的
[
k
,
i
−
2
]
[k, i-2]
[k,i−2] 并更新
[
k
,
j
]
[k, j]
[k,j] ,如下图所示。
1-9 行算法的复杂度为 O ( n 3 ) O(n^3) O(n3) 。于是,作者在每次循环 j j j 时加入了 BeamPrune( Q Q Q, j, b) ,伪代码如下所示
- function BeamPrune( Q Q Q, j, b)
- for [i, j] in Q do
- candidates[i] = Q 1 , i − 1 ⋅ Q i , j Q_{1, i-1} · Q_{i, j} Q1,i−1⋅Qi,j
- candidates = SelectTopB(candidates, b) # 选出前 b 个最大的值
- for [ i , j ] [i, j] [i,j] in Q Q Q do
- if i i i not in candidates then
- delete Q i , j Q_{i, j} Qi,j
这样就使得复杂度降低为
O
(
n
b
2
)
O(nb^2)
O(nb2) ,
b
b
b 通常取 100 。
核苷酸
i
i
i 与 核苷酸
j
j
j 配对的概率为
p
i
,
j
=
∑
y
∈
Y
(
x
)
,
(
i
,
j
)
∈
p
a
i
r
s
(
y
)
p
(
y
)
p
(
y
)
=
e
−
Δ
G
∘
(
y
)
R
T
Q
(
x
)
p_{i, j} = \sum_{y ∈ \mathcal{Y}(\textbf{x}), (i, j) ∈ pairs(y)} p(y) \\ p(y) = \frac{e^{-\frac{ΔG^{\circ}(y)}{RT}}}{Q(\textbf{x})}
pi,j=y∈Y(x),(i,j)∈pairs(y)∑p(y)p(y)=Q(x)e−RTΔG∘(y)
Limitations
- 本方法无法保证近似质量,因为 b 为固定值。但如果以阈值作为筛选条件的话又不能保证运行时间。
- LinearPartition 在随机序列,特别是存在很多冲突的配对选择时表现不好。