Constraint generation(CG) approach

Efficient parameter estimation for RNA secondary structure prediction

Year: 2007
Authors: Mirela Andronescu, Anne Condon, Holger H. Hoos, David H. Mathews, and Kevin P. Murphy
Journal Name: Bioinformatics

Motivation

1. 基于自由能的RNA二级序列预测模型Turner99，没有有效、快速的参数估计方法（使用大量结构数据和自由能数据训练）。
2. 数据集中的最小自由能(MFE)也许有噪声。

Research Objective

提出一种能应用于大样本的自由能参数计算方法（CG方法），且对自由能噪声具有鲁棒性。

Dataset

$(x,y_x)\in \mathcal{S}$ ，其中 $y_x$ 为序列 $x$ 的真实MFE结构。
$(x,y_x,e_x)\in \mathcal{T}$ ，其中 $e_x$ 为结构 $y_x$ 的自由能。

Background

已知序列 $x$ ，目标为预测具有最小自由能的结构 $y$ 。序列 $x$ 和结构 $y$ 的自由能计算公式为
$$\Delta G(x,y,\mathbf{\theta })=c(x,y{)}^T\mathbf{\theta }=\sum _{k=1}^K{c}_k(x,y){\theta }_k$$
其中 $K$ 是特征的个数， $c_k(x,y)$ 为特征 $k$ 在出现在结构 $y$ 中的次数， ${\theta }_k$ 是特征 $k$ 的能量参数。

CONTRAfold模型使用最大可能性(ML)方法来估计参数。定义在序列 $x$ 和自由能参数 $\mathbf{\theta }$ 确定的条件下，结构为 $y$ 的概率为
$$p(y|x,\mathbf{\theta })=\frac{1}{Z(x,\mathbf{\theta })}exp(-\frac{1}{RT}\Delta G(x,y,\mathbf{\theta }))$$
其中， $R$ 是气体常数， $T$ 是绝对温度， $Z(x,\mathbf{\theta })$ 为分段函数。
$p(y|x,\mathbf{\theta })$ 为凸函数，可以基于梯度优化概率对数 $L_{\mathcal{S}}(\mathbf{\theta })={\sum }_{(x,y_x)\in \mathcal{S}}logp(y_x|x,\mathbf{\theta })$ ，来估计参数 $\mathbf{\theta }$ 。
因为数据集中的 $e_x$ 可能存在噪声，作者在概率中加入了精度为 $\tau$ 的高斯分布，即
$$p(y|x,\mathbf{\theta })\propto L_{\mathcal{S}}(\mathbf{\theta })-\tau \sum _{(x,y_x,e_x)\in \mathcal{T}}(e_x-c(x,y_x{)}^T\mathbf{\theta }{)}^2$$
论文中原式为 $L_{\mathcal{S}}(\mathbf{\theta })+\tau {\sum }_{(x,y_x,e_x)\in \mathcal{T}}(e_x-c(x,y_x{)}^T\mathbf{\theta }{)}^2$ ，但高斯分布的指数项有 $-\frac{1}{2}$ ，所以我认为应把 $+$ 变为 $-$ 。直观上来看，当误差越小时概率越大。如理解有误请大佬指正。
这种方法存在两个问题

1. 算力要求过高
2. 该模型不能处理缺失（特征不完全）的数据集。

Method

非最优的结构 $y$ 的自由能大于MFE结构 $y_x$ 的自由能，约束表示为
$$\Delta G(x,y_x,\mathbf{\theta })<\Delta G(x,y,\mathbf{\theta })$$
其中， $(x,y_x)\in \mathcal{S}$ ， y ∈ Y x ∖ { y x } y ∈ Y_x \setminus \{y_x\} y∈Yx​∖{yx​} ， $Y_x$ 是序列 $x$ 的所有可能的二级结构集合。
因为数据集中的 $e_x$ 可能存在噪声，作者引入了松弛变量 ${\delta }_{x,y}\ge 0$ 。放松后的约束表示为
$$\begin{gathered} \Delta G(x,y_x,\mathbf{\theta })<\Delta G(x,y,\mathbf{\theta })+{\delta }_{x,y} \\ (c(x,y_x)-c(x,y){)}^T\mathbf{\theta }-{\delta }_{x,y}＜0 \\ {M}_{\mathcal{S}}\mathbf{\theta }-\mathbf{\delta }＜0 \end{gathered}$$
其中， $M_{\mathcal{S}}$ 为 $(c(x,y_x)-c(x,y){)}^T$ ， $\mathbf{\delta }$ 为 ${\delta }_{x,y}$ 向量。
优化问题表示为
$$\begin{array}{rl}& min\Vert \mathbf{\delta }{\Vert }_2^2\\ & s.t.\begin{array}{l}{M}_{\mathcal{S}}\mathbf{\theta }-\mathbf{\delta }＜0\\ \mathbf{\delta }\ge 0\end{array}\end{array}$$


作者又加入了数据集 $\mathcal{T}$ ，约束表示变为
$$\begin{gathered} \Delta G(x,y_x,\mathbf{\theta })-{\xi }_x=c(x,y_x{)}^T\mathbf{\theta }-{\xi }_x=e_x \\ {M}_{\mathcal{T}}\mathbf{\theta }-\mathbf{\xi }=\mathbf{e} \end{gathered}$$
其中， ${\xi }_x$ 为预测 $e_x$ 的误差， $M_{\mathcal{T}}$ 为 $c(x,y_x{)}^T$ ，$\mathbf{\xi }$ 为 ${\xi }_x$ 向量， $e$ 为 $e_x$ 向量。
优化问题改变为
$$\begin{array}{rl}& min(1-\lambda )\frac{1}{|\mathcal{S}|}\Vert {\mathbf{m}}^T\mathbf{\delta }{\Vert }_2^2+\lambda \frac{1}{|\mathcal{T}|}\Vert \mathbf{\xi }{\Vert }_2^2\\ & s.t.\begin{array}{l}{M}_{\mathcal{S}}\mathbf{\theta }-\mathbf{\delta }＜0\\ M_{\mathcal{T}}\mathbf{\theta }-\mathbf{\xi }=\mathbf{e}\\ \mathbf{\delta }\ge 0\end{array}\end{array}$$
其中， $|\mathcal{S}|$ 代表集合 $\mathcal{S}$ 的样本数， $m_x$ 为用于计算序列 $x$ 的 $M_{\mathcal{S}}$ 的约束个数的倒数。本人理解的约束个数为特征数量，即将 $\mathbf{\delta }$ 进行标准化，防止因特征数过多导致的 $\mathbf{\delta }$ 过大。 $\mathbf{m}$ 为 $m_x$ 的向量，超参数 $0\le \lambda \le 1$ 控制 $\mathcal{S}$ 和 $\mathcal{T}$ 的相对重要性。
如果某个特征很少甚至没有出现在数据集中，那么该特征所对应的能量参数 $\theta$ 就会失去限制。所以作者加入限制，即预测出的能量参数 $\theta$ 不能超出Turner99能量参数 ${\theta }_0$ 的一定范围，表示为
$${\theta }_0-B\le \theta \le {\theta }_0+B$$
其中， $B$ 为超参数。


上述优化问题的限制数量会随着输入尺寸指数性增加，因为集合 $Y_x$ 的尺寸随着数据集 $\mathcal{S}$ 的尺寸指数性增加。所以，作者采取迭代的方法，将 $M_{\mathcal{S}}\mathbf{\theta }-\mathbf{\delta }＜0$ 改为
$$(c(x,y_x)-c(x,y^{\prime }){)}^T{\theta }^{(i)}-{\delta }_{x,y^{\prime }}^{(i)}<0$$
其中， $y^{\prime }$ 为使用上一轮迭代参数 ${\theta }^{(i-1)}$ 预测的结构。预测方法用的是SimFold和Mfold。

Future Work

当特征覆盖率较小时，ML模型更加稳健。所以作者计划将CG模型和ML模型相结合，用ML估计少量不可靠参数，用CG估计剩余参数。作者还将探讨如何引入替代特征，例如同轴碱基对堆叠和多环未配对片段中的不对称性，来改进 RNA 二级结构预测。