CN114121146B - 一种基于并行和蒙特卡罗策略的rna三级结构预测方法 - Google Patents

Abstract

本发明公开了一种基于并行和蒙特卡罗策略的RNA三级结构预测方法，该方法属于结构预测领域。该方法包括利用并行机制进行构象空间采样；根据最新更新的能量函数进行打分；通过两轮势能评判对构象基于“Stepwise　ansatz”的蒙特卡洛操作进行合理性判断；最后经过结构完整性和建模精度的判断，对结果进行处理，直到获得稳定的高精度高完整性的RNA三级结构。本发明所提RNA三级结构预测方法能够获得高精度高完整性的RNA三级结构。基于并行和蒙特卡罗策略的RNA三级结构预测方法灵活性高，蒙特卡罗次数可指定，建模精度和建模时间成本可以由用户来衡量；该方法克服了现有技术方法对RNA模体建模不完整的问题；本发明增加了构象采样的广度和深度，减小了伪最低自由能的影响，提高建模精度。

Description

一种基于并行和蒙特卡罗策略的RNA三级结构预测方法

技术领域

本发明属于结构预测领域，具体涉及一种基于并行和蒙特卡罗策略的RNA三级结构预测方法。

背景技术

新的研究发现，RNA有一些复杂的生物学功能。结构决定功能，因此想要探究RNA的功能，就必须要先了解RNA的结构。目前国内外的RNA三级结构测定方法主要有两种。第一种方法是利用x射线、核磁共振以及冷冻电子电子显微镜等实验测定方法，采用实验的方法测得的结果比较精确且可靠，但是构象数量随着RNA长度的增加呈指数增长，导致成本较高。第二种方法是基于生物计算的结构预测方法，当前的RNA三级结构预测算法主要是基于知识挖掘的预测方法和基于物理的预测方法。基于知识挖掘的三级结构预测方法依赖已知的RNA模板库；基于物理的预测方法减少了对数据库的依赖，但是仍存在结构建模精度不够高的问题，无法满足当前的结构预测需要。因此针对这个现状，我们需要对现有方法进行革新。

在蛋白质领域，存在一个假设，假设大分子的原生构象有最低自由能，并且自由能函数近似为氢键、范德华力、静电力和溶剂化项之和。但是由于蛋白质和RNA分子折叠方式的不同，蛋白质的研究方法应用到RNA的研究中得到的结果不佳。因此我们针对现有技术的缺陷，依然假设大分子原生构象具有最低自由能，但为不同的三级相互作用分配不同的权值，线性加和得到自由能。此外，针对单线程构象能力受限制，我们采用了并行机制，同时对建模结果进行了多重判断，得到一个专门用于RNA三级结构预测的方法逐步蒙特卡罗并行化方法(SMCP)。

发明内容

针对现有RNA结构预测方法的缺陷，本发明提供了一种基于并行和蒙特卡罗策略的RNA三级结构预测方法SMCP。SMCP通过并行机制来增加构象采样的广度及深度，通过多重势能评判来对中间结果进行筛选，通过结果判断增加结果的完整性并提高建模精度。其目的在于解决当前RNA结构预测方法中单线程构象采样的缺陷以及建模完整度和精度不高的问题。

基于并行和蒙特卡罗策略的RNA三级结构预测方法，包括以下步骤：

(1)对RNA模体进行初始化，确定并行线程数n和蒙特卡罗次数m；n为大于1的自然数，m取200-50000的自然数；

(2)为RNA模体构建构象空间，利用并行机制和“Stepwise ansatz”假设对构象空间进行高效的构象采样，进行单个核苷酸的添加，删除，合并，重采样等操作，经过多次随机操作得到一个候选构象集；

(3)将步骤(2)得到的候选构象，利用能量函数进行势能值计算，生物分子势能值近似于Rosetta能量函数值，Rosetta能量值为根据公式△E_total＝∑_iω_iE_i(Θ_i,aa_i)计算所有能量项按权重缩放后的线性和，其中E_i是能量项，ω_i是每个能量项的权重，Θ_i是几何自由度，aa_i是化学恒等式；此外，计算过程中需要根据连接权重公式来计算各能量项的势能E_x，其中E_x是能量项x的势能值；

(4)对步骤(3)得到的构象势能值进行评判，经过步骤(2)的随机操作后，势能值的变化决定了对核苷酸的添加、删除，构象合并、重采样等操作是否可以被接受；根据标准：

决定对核苷酸的随机操作是否可接受，其中Metropolis标准由公式：表示，经过初步势能评判后得到真正的候选构象集；

(5)对步骤(4)得到的候选构象集进行进一步势能评判，具有低势能值的构象结构更稳定，因此综合所有线程选出势能值最低的构象作为当前最佳候选构象；

(6)对步骤(5)得到的当前最佳候选构象进行精度计算，RMSD是描述分子两种构象结构相似性的一个重要指标；根据公式来计算RMSD，以此来描述建模精度，其中/>是原子j与参考构象或m个等价原子的平均位置之间的距离；此外，通常会执行刚性叠加使RMSD最小化，然后返回最小值作为最终的精度值。根据公式

计算，其中n，ν表示给定的两个点；

(7)对步骤(5-6)得到的当前最佳候选构象进行精度判断；我们认为预测构象与实验测定的构象误差在以内时，预测构象就为原生构象(即要求建模精度/>)；故执行判断：

(8)对步骤(7)得到的当前最佳候选构象进行完整度判断，判断

(9)经过步骤(8)得到的构象是高精度高完整度的构象，得到最终的建模结果，随后利用UCSF Chimera进行可视化分析，可以通过UCSF Chimera对实验测定的构象和RNA三级结构预测方法预测到的构象进行对比分析。

优选地，步骤(1)中n的取值为3，m的取值为10000；

优选地，步骤(3)中所述获得构象势能值时包括以下步骤：

(1-1)计算原子对相互作用的能量。原子对间/内相互作用包括：范德华力，静电力，溶剂化项，氢键力，二硫键力。体现原子对相互作用的能量项包括：fa_rep，fa_intra_rep，fa_atr，fa_elec，fa_sol，lk_ball_wtd，hbond_sc，dslf_fal3，hbond_lr_bb，hbond_sr_bb，hbond_bb_sc；

(1-2)计算蛋白质骨架和侧链扭转相关的能量。表示扭转角的术语有：拉式图，骨架设计项，以及侧链构象，相关的能量项包括：rama_prepro，p_aa_pp，fa_dun；

(1-3)计算特殊情况下扭转项(肽键二面角)的能量。相关能量项包括:omega，pro_close，

yhh_plannarity；

(1-4)计算非理想键长和角度(笛卡尔积键能)的能量。相关的能量项包括：cart_bonded；

(1-5)Rosetta框架下所有能量函数的能量项均相同，不同能量函数的区别在于能量项权值的不同；

根据公式E_total＝ω_{fa_rep}E_{fa_rep}+ω_{fa_intra_rep}E_{fa_intra_rep}+ω_{fa_atr}E_{fa_atr}+ω_{fa_ elec}E_{fa_elec}E_{fa_elec}+ω_{fa_sol}E_{fa_sol}+ω_{lk_ball_wtd}E_{lk_ball_wtd}+ω_{hbond_sc}E_{hbond_sc}+ω_{dslf_fal3}E_{dslf_fal3}+ω_{hbond_lr_bb}E_{hbond_lr_bb}+ω_{hbond_sr_bb}E_{hbond_sr_bb}+ω_{hbond_bb_sc}E_{hbond_bb_sc}+ω_{rama_prepro}E_{rama_prepro}+ω_{p_aa_pp}E_{p_aa_pp}+ω_{fa_dun}E_{fa_dun}+ω_omegaE_omega+ω_{pro_close}E_{pro_close}+ω_{yhh_plannarity}E_{yhh_plammarity}+ω_{cart_bonded}E_{cart_bonded}计算上述步骤所有能量项的加权和，其中ωx是x能量项的权值，计算后得到候选构象的势能值。

优选地，步骤(4)核心在于计算系统能量变化ΔE。步骤(4)中所述Metropolis标准，根据公式来确定接受标准。其中df为新构象与原构象的适应度之差，即df＝f(new)-f(old)；t为退火过程的控制参数。

与现有技术方法相比，本发明的有益效果是：

本算法对RNA三级结构预测算法进行了革新，实现高效的结构预测。该算法基于“Stepwise ansatz”

假设，通过对单个核苷酸进行操作，避免了需要一次性枚举全部构象的情况；通过随机采样添加单个核苷酸后的构象来预测结构，实现了不依赖片段或粗粒度建模阶段，且减小了计算量，节省了建模时间；利用并行化对算法进行优化，多程序同时运行，层层根据能量值筛选，提高预测精度和建模完整度，节省建模时间。

附图说明

图1为并行机制的原理图；

图2为SMCP方法流程图；

图3为使用SMCP方法预测RNA三级结构的一个实例；

图4为使用SMCP方法与Rosetta框架下SWM方法预测RNA三级结构的建模精度结果对比图。

具体实施方式

为了清楚的阐明本发明的技术方案，下面结合附图(1-3)以及实例对本发明进行阐述，此处的实例仅用于解释本发明，并不限定本发明。

图1展示了串行采样和并行采样的原理图。采用串行采样方法时，从s开始随机搜索进行构象采样，通过蒙特卡罗机制可以找到局部最低能量的位置；然而，单线程的可搜索性是有限的，很难跨越能量障碍找到真正的最低能量，单线程构象搜索获得的构象的最低势能可能是伪最低势能，从而导致RNA三级结构预测方法的预测精度低。采用并行采样方法时，多个线程在不同的初始位置s开始随机搜索同一个构象空间，所有线程将获得一个局部最低能量谷；对所有线程采样获得的局部构象样本进行综合处理，增加了获得构象空间中实际最低能量谷的概率，从而获得高质量的样本提高预测精度。

图2展示了SMCP方法预测RNA三级结构的流程步骤。选定RNA模体实例为l1_sam_ll_riboswitch(PDB编号：2QWY，模体长度：7，序列：GCAGUCG)。SMCP方法的输入有两个pdb格式的3D结构文件，一个是l1_sam_ll_riboswitch模体的初始构象，SMCP方法将在此结构的基础上进行建模；另一个是l1_sam_ll_riboswitch模体的原生构象，原生结构即实验测定的结构，与SMCP方法预测得到的结构进行对比，用于结构预测方法的预测精度分析。此外，还需输入1个fasta序列文件，1个flags命令操作文件，并指定线程数n＝3，蒙特卡罗次数m＝10000。SMCP方法的输出为该方法预测的RNA三级结构及其结构预测精度。以下是RNA三级结构预测的具体步骤：

1、构象采样

使用并行机制和“Stepwise ansatz”假设对构象空间进行高效的构象采样(构象空间内包含7个核苷酸：GCAGUCG)，在已知GCUCG结构的前提下，进行核苷酸A和G的添加，删除，重采样等随机操作，经过10000次随机蒙特卡罗操作得到一个候选构象集，采样过程如下(仅举例)。

重采样成功操作(以第9904次采样为例)：

(1)建模:1-2 4-7

(2)重建模移动的与5号核苷酸U相连的4号核苷酸G

(3)RMSD 1.512(4号核苷酸G的23号原子),叠加到1-2 5-7号核苷酸的第86号原子上(RMSD0.0000007)

(4)尝试的次数：10000，成功的次数:13；

重采样失败操作(以第9999次采样为例)：

(1)建模:1-3 5-7

(2)重建模移动的与2号核苷酸C相连的3号核苷酸A

(3)RMSD 3.536(3号核苷酸A的22号原子),叠加到1-2 5-7号核苷酸的第86号原子上(RMSD0.0000005)

(4)尝试的次数：3092，成功的次数:20；

添加失败操作(以第9998次采样为例)，标准与重采样一致，故仅举失败的例子：

(1)建模:1-3 5-7；

(2)添加4号核苷酸G时，连接到3号核苷酸A上；

(3)RMSD 5.777(3号核苷酸A的27号原子),叠加到其它核苷酸的86号原子(RMSD0.0000008)；

(4)尝试添加位置的数：100000，成功的次数:17；

删除失败操作(以第10000次采样为例)，标准与重采样一致，故仅举失败的例子：

(1)建模:1-3 5-7；

(2)删除连接到2号核苷酸C上的3号核苷酸A；

(3)RMSD 0.000,叠加到其它核苷酸的86号原子(RMSD 0.0000003)；

2、能量函数打分

对采样得到的候选构象，利用能量函数进行势能值计算，生物分子势能值近似于Rosetta能量函数值，Rosetta能量值为根据公式△E_total＝∑_iω_iE_i(Θ_i,aa_i)计算所有能量项按权重缩放后的线性和，其中E_i是能量项，ω_i是每个能量项的权重，Θ_i是几何自由度，aa_i是化学恒等式；此外，计算过程中需要根据连接权重公式来计算各能量项的势能E_x，其中E_x是x能量项的势能值，计算过程包括以下步骤：

(1-1)计算原子对相互作用的能量。原子对间/内相互作用包括：范德华力，静电力，溶剂化项，氢键力，二硫键力。体现原子对相互作用的能量项包括：fa_rep，fa_intra_rep，fa_atr，fa_elec，fa_sol，lk_ball_wtd，hbond_sc，dslf_fal3，hbond_lr_bb，hbond_sr_bb，hbond_bb_sc；

(1-2)计算蛋白质骨架和侧链扭转相关的能量。表示扭转角的术语有：拉式图，骨架设计项，以及侧链构象，相关的能量项包括：rama_prepro，p_aa_pp，fa_dun；

(1-3)计算特殊情况下扭转项(肽键二面角)的能量。相关能量项包括:omega，pro_close，yhh_plannarity；

(1-4)计算非理想键长和角度(笛卡尔积键能)的能量。相关的能量项包括：cart_bonded；

(1-5)Rosetta框架下所有能量函数的能量项均相同，不同能量函数的区别在于能量项权值的不同；计算上述步骤所有能量项的加权和，得到候选构象的势能值。

经过能量函数的计算后，不同随机操作的势能变化如下：

重采样成功操作势能值变化(以第9904次采样为例)：-5.247→-7.460，势能值降低(初始结构势能值:-5.247)；

重采样失败操作势能值变化(以第9999次采样为例)：-2.184→-3.170，势能值降低(初始结构势能值:-2.184)；

添加失败操作势能值变化(以第9998次采样为例)：17.900→-1.671，势能值降低(初始结构势能值:15.326)；

删除失败操作势能值变化(以第10000次采样为例)：-3.702→-3.702，势能值不变(初始结构势能值:-3.702)；

3、势能评判进一步确定构象

势能评判的核心在于计算系统能量变化ΔE。势能值的变化决定了核苷酸的添加、删除，构象合并、重采样等操作是否可以被接受；

其中所述Metropolis标准，利用了蒙特卡罗的思想。能量升高时，生成0到1之间的一个随机数α，并与exp(-△E/kT)进行比较，如果α>exp(-△E/kT)则拒绝接受；反之，则接受，得到真正的候选构象集。

对每一个操作都要进行评判，最终选择接受或拒绝，势能评判过程如下：

重采样成功操作(以第9904次采样为例)：

(1)执行连接到5号核苷酸U的4号核苷酸G重采样的反操作：重采样连接到4号核苷酸G上的1-2 5-7号核苷酸；

(2)执行完成后,建模:1-2 4-7；

(3)势能值变化：-6.82358→-7.46016，势能值降低；

(4)蒙特卡罗操作是否被接受？接受(原操作和反操作势能值均降低)；

重采样失败操作(以第9999次采样为例)：

(1)执行连接到2号核苷酸C的3号核苷酸A重采样的反操作：重采样连接到2号核苷酸C上的3号核苷酸A；

(2)执行完成后,建模:1-3 5-7；

(3)势能值变化：-6.33765→-3.16991，势能值升高；

(4)蒙特卡罗操作是否被接受？拒绝(原操作势能值降低，反操作势能值升高)；

添加失败操作(以第9998次采样为例)：

(1)执行反操作，删除连接到3号核苷酸A上的4号核苷酸G；

(2)删除后,建模:1-7；

(3)势能值变化：-6.33765→-1.67202，势能值升高；

(4)蒙特卡罗操作是否被接受？拒绝(原操作势能值降低，反操作势能值升高)；

删除失败操作(以第10000次采样为例)：

(1)执行反操作，添加连接到2号核苷酸C上的3号核苷酸A；

(2)执行完成后,建模:1-2 5-7；

(3)势能值变化：-6.33765→-3.70202，势能值降低(初始结构势能值:-3.702，与初始值无变化)；

(4)蒙特卡罗操作是否被接受？拒绝(原操作势能值不变，反操作势能值降低，但与初始值一致)；

4、多线程综合势能评判

具有低势能值的构象结构更稳定，因此综合所有线程选出势能值最低的构象作为当前最佳候选构象。针对l1_sam_ll_riboswitch模体所做的3个线程所得的构象势能值分别为：-11.133REU(REU：Rosetta Energy Units)，-10.123REU，-12.155REU，根据原则：结构势能值越低，结构越稳定，选择势能值最低的构象，即势能值为-12.155REU的构象。

5、建模精度计算

RMSD是描述分子两种构象结构相似性的一个重要指标；根据基本公式来计算RMSD，以此来描述建模精度，其中δ_i是原子i与参考构象或n个等价原子的平均位置之间的距离；计算过程中通常会执行一个刚性叠加，使RMSD最小化，然后返回这个最小值作为最终的精度值RMSD。此时需要根据公式来计算RMSD，其中n，ν表示给定的两个点。通过上述计算公式可得l1_sam_ll_riboswitch模体建模精度为/>其势能值为-12.155REU。

6、建模精度和完整度判断及处理

对当前最佳候选构象进行精度判断；判断当前最佳构象建模精度为可满足该精度要求，无需重新建模；

对当前最佳候选构象进行完整度判断，判断

当前构象的missing值为0，表明SMCP方法已经完成了对GCAGUCG这7个核苷酸的完整建模；

对l1_sam_ll_riboswitch模体进行10000次建模完成后，蒙特卡罗随机操作的相关统计如下：

(1)添加次数：1095；接受率：0.2868；

(2)删除次数：3968；接受率：0.0769；

(3)重采样次数：4937；接受率：0.4588；

经过随机蒙特卡罗并行化采样，以及对精度和完整度判断和处理，将A和G核苷酸添加到了A链上的第24和第25号位置上，该构象即为最终的建模结果。

7、构象可视化分析

建模后得到的高精度高完整度的结构，利用UCSF Chimera进行可视化分析，可以通过UCSF Chimera对实验测定的构象和SMCP方法预测到的构象进行对比分析。对比结果如图3所示，其中图A为l1_sam_ll_riboswitch模体的实验测定结构，图B为SMCP方法预测的l1_sam_ll_riboswitch模体的结构，从图3中A与B图可以看出，SMCP方法预测的RNA结构与实验测定的真实结构相似度极高；从上述建模结果数据来看，SMCP方法对l1_sam_ll_riboswitch模体进行建模的RMSD为而利用Rosetta框架下RNA三级结构预测当前最佳方法SWM得到的RMSD为/>表明SMCP方法比SWM方法预测RNA三级结构精度高。

利用SMCP和Rosetta框架下SWM方法对由9个RNA组成的基准进行建模，图4为使用SMCP方法与SWM方法对上述基准进行三级结构建模的RMSD结果对比图，图中横坐标为使用SWM方法预测RNA三级结构的RMSD，纵坐标为使用SMCP方法预测RNA三级结构的RMSD，图中每一个点/方形代表一个RNA模体。

从图4可以看出，对基准中的9个RNA模体进行结构预测时，SMCP方法建模得到的RMSD值更低，即建模精度更高，这表明SMCP方法在RNA三级结构高精度建模领域占有优势地位，SMCP方法预测RNA三级结构时预测精度更高；

图4中存在2个黑色方形标识的RNA模体，利用SWM方法进行RNA三级结构预测时，无法实现这对2个RNA模体的完整建模；而利用SMCP方法可以实现对这2个RNA模体中所有核苷酸的完整建模，这表明SMCP方法预测RNA三级结构时预测完整度更高。

Claims (2)

1.一种基于并行和蒙特卡罗策略的RNA三级结构预测方法，其特征在于，包括如下步骤：

（1）对RNA模体进行初始化，确定并行线程数n和蒙特卡罗次数m；n为大于1的自然数，m取200-50000的自然数；

（2）为RNA模体构建构象空间，利用并行机制和“Stepwise ansatz”假设对构象空间进行高效的构象采样，进行单个核苷酸的添加，删除，合并，重采样操作，经过多次随机操作得到一个候选构象集；

（3）将步骤（2）得到的候选构象，利用能量函数进行势能值计算，生物分子势能值近似于Rosetta能量函数值，Rosetta能量值为根据公式计算所有能量项按权重缩放后的线性和，其中/>是能量项，/>是每个能量项的权重，/>是几何自由度，是化学恒等式；此外，计算过程中需要根据连接权重公式/> 来计算各能量项的势能；

所述步骤（3）中所述获得构象势能值时包括以下步骤：

（1-1）计算原子对相互作用的能量：原子对间/内相互作用包括：范德华力，静电力，溶剂化项，氢键力，二硫键力；体现原子对相互作用的能量项包括：，/>，，/>，/>，/>，/>，/>，/>，，/>；

（1-2）计算蛋白质骨架和侧链扭转相关的能量：表示扭转角的术语有：拉式图，骨架设计项，以及侧链构象，相关的能量项包括：，/>，/>；

（1-3）计算扭转项肽键二面角的能量：相关能量项包括：，/>，；

（1-4）计算非理想键长和角度笛卡尔积键能的能量：相关的能量项包括：；

（1-5）Rosetta框架下所有能量函数的能量项均相同，不同能量函数的区别在于能量项权值的不同；根据公式

计算上述步骤所有能量项的加权和；

（4）对步骤（3）得到的构象势能值进行评判，经过步骤（2）的随机操作后，势能值的变化决定了核苷酸的添加、删除，构象合并、重采样操作是否可以被接受；根据标准：，决定对核苷酸的随机操作是否可接受，其中Metropolis标准由公式：表示，经过初步势能评判后得到真正的候选构象集；α为一个0-1之间的随机数；

（5）对步骤（4）得到的候选构象集进行进一步势能评判，具有低势能值的构象结构更稳定，因此综合所有线程选出势能值最低的构象作为当前最佳候选构象；

（6）对步骤（5）得到的当前最佳候选构象进行精度计算，RMSD是描述分子两种构象结构相似性的一个重要指标；根据公式来计算RMSD，以此来描述建模精度，其中/>是原子j与m个等价原子的平均位置之间的距离；此外，会执行刚性叠加使RMSD最小化，然后返回最小值作为最终的精度值；根据公式计算，其中/>，/>表示给定的两个点；

（7）对步骤（5）-（6）得到的当前最佳候选构象进行精度判断；预测构象与实验测定的构象误差在2Å以内时，预测构象就为原生构象，即要求建模精度RMSD2Å；

故执行判断：；

（8）对步骤（7）得到的当前最佳候选构象进行完整度判断，判断；

（9）经过步骤（8）得到的构象是高精度高完整度的构象，得到预测的最终的RNA三级结构。

2.根据权利要求1所述的一种基于并行和蒙特卡罗策略的RNA三级结构预测方法，其特征在于，步骤（1）中n的取值为3，m的取值为10000。