2. 苏州大学江苏省计算机信息处理技术重点实验室,苏州 215006
2. Provincial Key Laboratory for Computer Information Processing Technology of Jiangsu, Suzhou 215006, China
近年来,从头预测的方法在蛋白质结构预测领域取得了不错的成绩。但是,就目前的技术手段而言,从头预测的方法仍然面临两大基本难点[1]:第一,由于蛋白质内部原子相互作用的复杂性,人们难以找到足够准确的能量函数来描述一个蛋白质构象;第二,蛋白质结构的构象空间相当大,尤其在残基序列较长的时候,如果没有合适的、高效的采样算法,蛋白质构象采样将是一个“灾难性”的计算问题。由此可见,采样算法在从头预测的方法中有着举足轻重的地位[2],因此,许多研究者针对这一问题开发了多种适用于蛋白质三维结构采样的算法,比如改进的遗传算法[3]、构象空间退火[4-5]和分子动力学模拟采样[6]。
在蛋白质结构预测的语境下,采样是指产生一个满足要求的构象,而根据这个要求的不同,我们可以将蛋白质构象采样分为狭义的采样和广义的采样:
(1) 狭义的采样:一个蛋白质构象可以通过不同的自由度来表示,常见的有二面角、原子坐标等。如果使用二面角表示蛋白质构象,那么这里的采样就是指对二面角的一次赋值;如果使用原子坐标来表示蛋白质构象,那么这里的采样就是指对原子坐标的一次赋值。Rosetta[7]中的每一个Mover都是这种意义下的采样。
(2) 广义的采样:对于表示蛋白质构象的自由度,如果我们认为它们满足某一特定的概率分布,比如正态分布N(d, u),那么就可以对这个特定的概率分布进行一次采样,同样可以得到该自由度表示下的蛋白质的构象,而这样的采样就是广义上的采样。
对于狭义的采样,由于是在自由度上进行不断的赋值尝试,再通过能量函数进行打分评估。不可避免的问题就是随着残基数量的增加,自由度的数量会成指数级增长,即使在现有的计算能力下,这样的计算量增长对采样算法来说,仍然是不可行的。而且,能量函数本身就是一个值得研究的问题。本文提出的基于距离约束模型的HMC(Hybrid Monte Carlo)采样算法则是对概率模型的采样,即使残基数量增加,蛋白质序列变长,但是在参数数量上面临的问题会比直接对自由度赋值的方法好许多。
HMC采样[8-9]是一种对概率模型的采样算法,其采样机制通过引入Hamiltonian动力学系统[10],避免在巨大的采样空间中随机采样,并且可以在采样的过程中保留概率分布所具有的统计规律。例如,在本文的案例中,如果我们使用原子坐标作为自由度来表示一个蛋白质结构,每个残基选取4个原子,即N、Cα、C和O;每个原子有3个笛卡尔坐标,即x、y和z坐标。因此,表示一个残基的就是12个自由度。对于一个长度为10残基的蛋白质构象,用来表示它的自由度就有120个。进一步认为这120个自由度满足一个120个随机变量的概率分布,那么就可以用HMC算法对这个概率分布进行一次采样,可以认为这是对该结构的一次“理性”的采样。因为采样的过程是模拟的Hamiltonian动力学系统,即粒子在空间中运行的状态,那么如何保证蛋白质结构在采样的过程中不会因为随机的走动而散架呢?本文引入了距离约束模型。
蛋白质的空间结构并不仅仅是由一个个的残基的独立位置决定的。事实上,残基间的相互作用对蛋白质结构有着巨大的影响[11]。在蛋白质结构中,按照残基的间隔,一般可以分为三类距离约束[12]: 1) 短距离约束(Short-range),此类约束的残基间隔相对较短,范围通常在6-12个残基;2) 中等距离约束(Medium-range),此类约束的残基间隔一般在12-24个残基;3) 长距离约束(Long-range),此类约束的残基间隔一般在24个残基以上。同时,在蛋白质结构中,一般认为残基空间欧氏距离在8 Å以内的残基间是有相互作用的,因此本次实验中使用的约束模型就是8 Å以内的长距离约束。
1 材料与方法 1.1 数据准备蛋白质结构预测评估大会[13](CASP: Critical Assessment of techniques for protein Structure Prediction)是世界范围内对蛋白质结构预测技术进行评估,评价当前蛋白质结构建模方法的能力及局限性的竞赛,被称为蛋白质结构预测领域的“奥运会”。本文从CASP8、CASP9和CASP10当中共选取了16个案例,它们的平均残基长度为100~200,最短的包含101个残基,最长的包含195个残基,每个蛋白质都是提取的单链。
如上文所述,本文采用的自由度是原子三维坐标,每个残基取其包括骨架原子Cα在内的4个原子的三维坐标,即每个残基由12个值表示。每条蛋白质为一行,以300行为一个batch。因为本文的主要目的是验证采样算法的可行性,所以在数据准备的时候是从蛋白质的天然结构里获取其距离约束的,而在实际应用的时候可以考虑使用同源蛋白质的距离约束。
从天然结构中获取距离约束的过程如下:首先,遍历天然结构的残基,计算出每对残基间隔大于24个残基的残基对的空间欧氏距离,然后找出小于8 Å的残基对,记录下这些残基对以及距离值,作为约束。采样的输入数据是用Rosetta平台ab initio得到的。考虑到ab initio得到的结构具有较大的RMSD,可以更加显著的看到采样算法的效果,而且,采样是为了得到更加符合距离约束的decoys,因此选用这样的数据集作为采样的输入。
1.2 HMC采样HMC采样算法的核心在于通过引入Hamiltonian动力学系统,模拟粒子在空间的运动,从而对一个能量模型进行采样。事实上,蛋白质的空间结构就是由空间中一个个粒子组成的,对蛋白质结构的采样可以想象为让这些粒子在空间中不断的运动,以找到和天然结构最接近的构象。可见,如果不添加任何约束,那么采样空间将接近于无穷大,因为粒子在空间里可以向任意方向以任意速度运动。面对如此巨大的采样空间,任何采样方法都是无能为力的。因此必须添加一定的约束从而一定程度上缩小采样空间。
1.2.1 Hamiltonian动力学系统粒子在空间中的任何时刻都具有两种属性,即速度和位置,这两种属性对应于两种能量:动能和势能,我们把这样的一个动力系统称作Hamiltonian动力学系统。HMC采样正是通过模拟这样的一个系统,对一个高维空间进行采样的。其能量函数的定义如公式(1)、(2)、(3)所示:
$E\left( s \right) = s \in {R^D}$ | (1) |
$K\left( \varphi \right) = \varphi \in {R^D}$ | (2) |
$H\left( {s, \varphi } \right) = E\left( s \right) + K\left( \varphi \right) = E\left( s \right) + \frac{1}{2}\sum\limits_i {\varphi _i^2} $ | (3) |
其中,E(s)为势能,K(φ)为动能。在本文的案例中,E(s)即为距离约束模型的值,在采样过程中,会沿着这个能量和的最小方向进行采样。这样的过程就一定程度上缩小了采样空间,避免了直接采样的不可行性。而在Hamiltonian系统中,我们并没有直接对p(s)进行采样,所用的方法是通过对一个正则分布:
$p\left( {s, \varphi } \right) = \frac{1}{Z}\exp \left( {-H\left( {s, \varphi } \right)} \right) = p\left( s \right)p\left( \varphi \right)$ |
进行采样,这是因为两个概率分布p(s)和p(φ)是相互独立的,边缘化φ并不会产生很大的影响,并且可以重现原有的概率分布。
1.2.2 蛙跳算法因为模拟的时间是离散的,而粒子的运动是连续的,所以在实现过程中并不能很精确的模拟这样一个系统。对于这个问题,有多种方案可以解决,其中比较好的是蛙跳算法[14-15],可以通过三步操作解决这个问题,如公式(4)、(5)、(6)所示:
${\varphi _i}\left( {t + \frac{\varepsilon }{2}} \right) = {\varphi _i}\left( t \right)-\frac{1}{2}\frac{\partial }{{{\partial _{{s_i}}}}}E\left( {s\left( t \right)} \right)$ | (4) |
${s_i}\left( {t + \varepsilon } \right) = {s_i}\left( t \right) + \varepsilon {\varphi _i}\left( {t + \frac{\varepsilon }{2}} \right)$ | (5) |
${\varphi _i}\left( {t + \varepsilon } \right) = {\varphi _i}\left( {t + \frac{\varepsilon }{2}} \right)-\frac{\varepsilon }{2}\frac{\partial }{{{\partial _{{s_i}}}}}E\left( {s\left( {t + \varepsilon } \right)} \right)$ | (6) |
这里首先计算了速度在时间为t+ε/2时候的势能,然后通过这个值计算出si(t+ε)和φi(t+ε)。
1.2.3 Metropolis接受/拒绝由于使用有限的步长ε并不能准确的表达能量H(s, φ),因此可能在模拟的过程中产生误差。而且使用单精度的数值也可能导致舍入误差,所以对于一次新的采样,会使用Metropolis准则[16-17]进行一次接受/拒绝判定。其判定过程如公式(7)所示:
${P_{acc}}\left( {x, x'} \right) = \min \left( {1, \frac{{\exp \left( {-H\left( {s', \varphi '} \right)} \right)}}{{\exp \left( {-H\left( {s, \varphi } \right)} \right)}}} \right)$ | (7) |
Metropolis准则是模拟退火算法中最常用的接受准则,其基本思想是以一定的概率接受新状态。它的重点在于如果H(s′, φ′)小于H(s, φ)那么,新状态被接受,反之,则以Pacc(x, x′)的概率接受,而不是拒绝。
1.2.4 HMC采样流程在本文的实验中,势能是根据蛋白质残基的距离约束构建的模型,动能是一个随机的均匀分布,整个采样的流程如下:
(1) 从均匀分布中产生一个随机的速度(动能);
(2) 执行n次蛙跳算法,获取新的状态x′(势能+动能),即对残基的三维坐标进行n次调整,模拟粒子的运动,每次赋值都是通过计算势能函数的梯度得到的,并不是随机的赋值;
(3) 对新的状态x′和旧的状态χ执行一次Metropolis接受判断,是否可以接受;
(4) 根据当前接受率和目标接受率对步长进行一次调整,调整采样效率;
(5) 重复以上步骤,直到达到特定的步数;
HMC采样的过程中,有几个重要的参数,比如步长、蛙跳次数、步长增益幅度等等,这些参数直接涉及到采样的效果和速度(见表 1)。
实验中,步长的初始值和变化范围设置的比较小,这是可以理解的,因为蛋白质原子空间的变化本身就不是很大,太大幅度的变化可能导致空间严重变形,或者直接越过最佳位置。而且本文仅仅加了距离约束,没有添加其它约束,如果设置过大,那么将会出现不可预期的结构。
1.3模型构建和实验流程
实验中使用的是残基间隔大于24并且残基空间欧氏距离小于8 Å的残基对作为约束。在实验的过程中,对于每一次的采样结果,都会与天然结构进行计算。主要分为三步:
(1) 计算标记的每个残基对的空间欧氏距离;
(2) 将这些距离与天然结构中对应残基对的空间欧氏距离作差;
(3) 计算这些差值的最小平方和。
用这个值作为HMC采样的势能函数。这里,使用最小平方和作为模型输出,是因为该方法本身具有的优点:它可以找出与数据最佳的匹配函数,并且对于曲线有很好的拟合能力。在本实验中,该模型如公式(8)~(11)所示:
${\Gamma _D} = \sqrt {\left( {{{\left( {{D_{lx}}-{D_{rx}}} \right)}^2} + {{\left( {{D_{ly}} + {D_{ry}}} \right)}^2} + {{\left( {{D_{lz}}-{D_{rz}}} \right)}^2}} \right)} $ | (8) |
${\Gamma _N} = \sqrt {\left( {{{\left( {{N_{lx}}-{N_{rx}}} \right)}^2} + {{\left( {{N_{ly}} + {N_{ry}}} \right)}^2} + {{\left( {{N_{lz}}-{N_{rz}}} \right)}^2}} \right)} $ | (9) |
其中,ΓD和ΓN分别表示decoys和天然结构中,两个标记的残基对之间的空间欧氏距离。然后对它们求差值:
${\Delta _i} = {\Gamma _{Di}}-{\Gamma _{Ni}}$ | (10) |
最后,再将这些差值求平方和:
$\varphi = \sum\limits_{i = 0}^n {{\Delta _i}} $ | (11) |
φ的大小反映了采样得到的构象是否满足从天然结构中获取的距离约束。其值越小,说明标记的每对残基对的空间欧氏距离和天然结构中对应的残基对的空间欧氏距离越接近,因此我们也就认为在采样的过程中,距离约束起到了作用,得到了更加满足距离约束的构象。
按照上述模型,在本地机器上(CPU:Intel(R) Xeon(R) E5-2620 v2 @ 2.10Hz内存:16GB)跑完HMC算法用时3.5 h,每个案例在6~26 min不等。整个采样实验流程如下:
(1) 使用Rosetta平台,对每个案例用ab initio的方法得到300个decoys作为采样的输入;
(2) 从天然结构中计算得到距离约束,即包括残基对和对应的空间欧氏距离;
(3) 对于ab initio得到的decoys,提取每个残基骨架原子的三维坐标,构成一个行数为300,列数为序列长度×12的输入矩阵;
(4) 将做好的矩阵输入到HMC当中,进行多次的迭代采样,并输出最终得到的构象。
之所以选用Rosetta平台制作输入集,也是为了作为对比结果。Rosetta的ab initio就是一种不带距离约束的采样,本文使用这样的decoys作为输入,经过带距离约束的HMC采样,然后比较输出decoys的距离约束,分析HMC采样的效果。
2 结果与分析实验一共使用了16个案例,每个案例用300个decoys作为采样的输入,每个案例采样1 001次,其中有1 000次作为HMC的burn-in阶段,丢弃。保留最后一次的结果。对比采样之前的input-decoys和采样之后得到的output-decoys,分别对每个decoys中标记的距离小于8Å的残基对的个数变化做了统计并比较,变化比例的计算如公式(12)所示:
$\begin{array}{c} 采样前大于8Å\\ 变化比例 = \frac{{采样后小于8Å的残基数}}{{总距离约束数量}} \end{array}$ | (12) |
实验统计结果如表 2所示,其中第一列是所选案例的PDB ID,二、三两列分别是残基和约束对的个数,最后三列则是变化的百分比。
从表 2中可以看出,在所有的案例中,采样后的300个decoys里,小于8Å的残基对个数有了明显的增长。如案例3JV6里,最大的一个构象有75.56%的残基对在采样之后空间欧氏距离小于8Å。即3JV6的天然结构中,一共有104个残基,长距离残基对个数为:1+2+3+…+79+80=3 240。其中,空间欧氏距离小于8Å的有132个。变化最大的一个decoy,在采样的过程中,这132个残基对里有132×0.755 6≈100个残基对,由原来大于8Å的距离变成小于8Å的距离,由不满足距离约束变为了满足距离约束。这就可以说明,本次实验中,HMC采样成功的加入了距离约束,使一些原本大于8Å的残基对的距离小于8Å。这样的意义在于使原来没有满足距离约束的残基对变得满足距离约束,如此一来,这样的构象就是比较好的构象,对下一步的三级结构预测有很大的意义。
另外,表 2是根据残基数从少到多的顺序排列的,通过观察数据还可以看出,残基数对采样效果的影响很少。残基数较少的案例其采样的效果并不一定好,例如1WWV、1WWW等;而残基数较多的案例其采样的效果也不一定会差,例如1HH4、2G0N、2WKP。这也印证了之前所说的,HMC采样避免了因序列增长,自由度数量急剧增长而导致的采样空间巨大,采样算法不可行的问题。
同时,本文也计算了采样前后蛋白质结构RMSD的变化,发现RMSD的变化很小,分析可能的原因有两点:
第一,采样的步长取值很小,也就意味着每一次采样对三维坐标的变化很小。模拟到Hamiltonian动力学系统中就是指粒子在空间中每次移动的速度很慢,因此原子移动距离很小,而且,实验中仅仅取了Cα,并没有全原子采样。事实上,如果步长较大,那么反而会出现不可预期的结果。虽然RMSD变化不大,但是在经过多次移动之后,仍然可以将大于8 Å的残基距离缩小到8 Å以内。可见采样的过程还是符合预期的结果;
第二,在采样的过程中,标记出的约束对里,虽然有部分残基对空间欧氏距离从大于8 Å降到8Å以内,但也有部分从8Å以内拉长到8Å以外。但是,从整体相互作用的残基对数变化的情况来看,采样还是可以获得比较好的构象的。事实上,采样得到的decoys并不是都要的,我们只需要保留比较好的构象。本文对这部分数据也做了统计,如表格3所示。
从表 3中可以看出,对于好的采样案例,空间欧氏距离小于8Å的残基对个数是增加的,而且增加的不少。在统计的16个案例中,有3个案例出现了超过50%的增幅(3JV5、3JV6、2WKP)。表中还统计了前1%平均和前10%平均,可以说,这些构象是比较好的构象,也正是我们需要的。结合表 2和表 3可见,基于距离约束的HMC采样算法在蛋白质结构预测方面是可行的。并且,相对于目前常用的Rosetta方法,能够得到更加满足距离约束的采样结果,取得更好效果,对于三级结构预测有很大的帮助。
本文以蛋白质结构中的距离约束为模型,使用HMC采样对蛋白质空间三维坐标进行采样。实验已可以证明HMC采样可以用在蛋白质结构预测方面,并且相对于目前常用的Rosetta有不错的效果。因为模拟了Hamiltonian动力学系统,避免了其它采样的随机性和盲目性,所以即使同时加入了上百对距离约束,仍然能在计算能力可及的范围内取得结果。然而,实验本身并不是完美的,仍然有需要改进的地方。比如,在加入距离约束的同时,如何能降低蛋白质结构的RMSD。因为结构预测的最终目的是获得近天然结构,距离约束只是为了更好的预测结构。因此,如何优化模型,加入其它信息,是本文的下一个工作重点。
[1] | HARDIN C, POGORELOV T V, Luthey-Schulten Z. Ab initio protein structure prediction[J]. Current Opinion in Structural Biology, 2002, 12(2): 176–181. DOI:10.1016/S0959-440X(02)00306-8 (0) |
[2] | KIM D E, BLUM B, BRADLEY P, et al. Sampling bottlenecks in de novo protein structure prediction[J]. Journal of Molecular Biology, 2009, 393(1): 249–260. DOI:10.1016/j.jmb.2009.07.063 (0) |
[3] | JUDSON R S, COLVIN M E, MEZA J C, et al. Do intelligent configuration search techniques outperform random search for large molecules?[J]. International Journal of Quantum Chemistry, 1992, 44(2): 277–290. DOI:10.1002/(ISSN)1097-461X (0) |
[4] | LEE J, SCHERAGA H A, RACKOVSKY S. New optimization method for conformational energy calculations on polypeptides: conformational space annealing[J]. Journal of Computational Chemistry, 1997, 18(9): 1222–1232. DOI:10.1002/(ISSN)1096-987X (0) |
[5] | LEE J, SCHERAGA H A. Conformational space annealing by parallel computations: extensive conformational search of Met-enkephalin and of the 20-residue membrane-bound portion of melittin[J]. International Journal of Quantum Chemistry, 1999, 75(3): 255–265. DOI:10.1002/(ISSN)1097-461X (0) |
[6] | LEE M R, TSAI J, BAKER D, et al. Molecular dynamics in the endgame of protein structure prediction[J]. Journal of Molecular Biology, 2001, 313(2): 417–430. DOI:10.1006/jmbi.2001.5032 (0) |
[7] | ROHL C A, STRAUSS C E M, MISURA K M S, et al. Protein structure prediction using Rosetta[J]. Methods in Enzymology, 2004, 383: 66–93. DOI:10.1016/S0076-6879(04)83004-0 (0) |
[8] | DUANE S, KENNEDY A D, PENDLETON B J, et al. Hybrid monte carlo[J]. Physics Letters B, 1987, 195(2): 216–222. DOI:10.1016/0370-2693(87)91197-X (0) |
[9] | SEXTON J C, WEINGARTEN D H. Hamiltonian evolution for the hybrid monte carlo algorithm[J]. Nuclear Physics B, 1992, 380(3): 665–677. DOI:10.1016/0550-3213(92)90263-B (0) |
[10] | DIRAC P A M. Generalized hamiltonian dynamics[C] //Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. The Royal Society, 1958, 246(1246): 326-332. (0) |
[11] | WANG Z, XU J. Predicting protein contact map using evolutionary and physical constraints by integer programming[J]. Bioinformatics, 2013, 29(13): i266–i273. DOI:10.1093/bioinformatics/btt211 (0) |
[12] | MONASTYRSKYY B, DANDREA D, FIDELIS K, et al. Evaluation of residue-residue contact prediction in CASP10[J]. Proteins: Structure, Function, and Bioinformatics, 2014, 82(S2): 138–153. (0) |
[13] | KRYSHTAFOVYCH A, MONASTYRSKYY B, FIDELIS K. CASP prediction center infrastructure and evaluation measures in CASP10 and CASP ROLL[J]. Proteins, 2014, 82(S2): 7–13. (0) |
[14] | VAN GUNSTEREN W F, BERENDSEN H J C. A leap-frog algorithm for stochastic dynamics[J]. Molecular Simulation, 1988, 1(3): 173–185. DOI:10.1080/08927028808080941 (0) |
[15] | SNYMAN J A. The LFOPC leap-frog algorithm for constrained optimization[J]. Computers & Mathematics with Applications, 2000, 40(8): 1085–1096. (0) |
[16] | NEAL R M, EDU E R T. Probabilistic inference using markov chain monte carlo methods[J]. Physics and Chemistry, Nato Science Series C: Mathematical and Physical Sciences, 1993, 92(440): 497–537. (0) |
[17] | ULYANOV N B, SCHMITZ U, JAMES T L. Metropolis monte carlo calculations of DNA structure using internal coordinates and NMR distance restraints:an alternative method for generating a high-resolution solution structure[J]. Journal of Biomolecular NMR, 1993, 3(5): 547–568. (0) |