一种统计计算固有无序蛋白多态结构的方法与流程

本发明涉及统计模拟与分析固有无序蛋白在折叠或其他非平衡动力学过程中形成的各种结构聚类、结构图景、结构特性等，尤其涉及一种统计计算固有无序蛋白多态结构的方法，属于蛋白质非平衡动力学与多态结构研究领域。

背景技术：

固有无序蛋白是一种缺乏固定或有序三维空间结构的蛋白质，它通常具有一系列、从完全无结构化到部分结构化(包括无规则卷曲、熔球态)的构象状态。传统的蛋白质结构与功能关系范例显示：蛋白质的功能依赖于固定的三维空间结构。然而，固有无序蛋白的许多研究结果对此提出了挑战，它们在生物体内依据不同的细胞环境会调整成不同的空间构型，形成一个结构或构象聚类。尽管缺乏稳定的结构，固有无序蛋白仍是一大类非常重要的功能蛋白，它在氨基酸序列、结构、功能、相互作用、动力学演化等方面具有不同于其他类型蛋白的特性。固有无序蛋白可形成蛋白质的几种典型结构类型，例如球状、纤维状和膜蛋白。固有无序蛋白的错误折叠、低聚和纤维状沉积是癌症、心血管疾病、神经退行性疾病等许多人类重大疾病的致病因素，研究其中间结构图景和特性对阐明这些疾病的致病机理具有重要的基础和应用意义。

实验上，X-Ray晶体衍射获取固有无序蛋白的电子密度谱还很缺乏，NMR是探测这类蛋白质结构和动力学信息的主要方法。然而，目前的NMR研究结果仅给出了固有无序蛋白的某个中间瞬态结构或与其他蛋白绑定后的构型等有限的信息，还没有关于这类蛋白质的结构图景或结构聚类的报道。理论上，蛋白质分子量大、结构复杂，计算时还需考虑蛋白质所处的环境因素，如溶剂、温度和压力等。采用传统的全原子或粗粒化分子动力学技术模拟蛋白质的折叠与非平衡动力学过程，计算量非常大。国内外已有的方法，如Gromacs、Amber、UNRES等，仅限于模拟很短、快速折叠的蛋白质的折叠过程。许多固有无序蛋白的折叠与动力学过程周期长、过程复杂，而且探查这类蛋白的结构图景或聚类需要进行大量的统计计算，分子动力学方法难以实现。到目前，还没有这方面的理论方法和结果。

技术实现要素：

本发明的目的是针对目前还没有固有无序蛋白结构图景或结构聚类的统计计算技术现状，提出了一种统计计算固有无序蛋白多态结构的方法。

本发明所提方法通过对固有无序蛋白进行上千次的折叠或非平衡动力学模拟，统计分析获得固有无序蛋白的结构图景、聚类和特性为主的重要信息。

为实现上述目的，一种统计计算固有无序蛋白多态结构的方法，步骤如下：

步骤(1)：获取蛋白质实验结构；

其中，蛋白质实验结构过程可以从蛋白质数据银行(PDB，http://www.rcsb.org)获取，且获取固有无序蛋白；蛋白质实验结构可以是X-Ray晶体衍射或NMR解析的结构；

步骤(2)：采用步骤(1)的蛋白质实验结构中心碳原子坐标数据，计算中心碳原子链接形成的键角和扭转角；

其中，中心碳原子记为Cα，中心碳原子的坐标记为(x,y,z)，键角与扭转角，记为(κ⁰i,τ⁰i)；

i从蛋白质的氮末端到碳末端依次编号，代表从氮末端开始的第i个氨基酸；

步骤(3)：构建蛋白质拟合模型；

其中，蛋白质等同于步骤(2)中的中心碳原子链；

步骤(3)具体为：

步骤(3).1根据步骤(2)输出将中心碳原子链分为多个片段，并对每个片段搜索能量函数的极小值，得到每个片段的拟合键角与扭转角及相应能量函数的参数值；

步骤(3).2将步骤(3).1中的每个片段顺次拼接在一起，再次搜索能量函数的全局极小，获得若干拟合结构及最优结构；

其中，步骤(3).1和步骤(3).2中的能量函数表述为如下公式(1)：

对上述公式(1)取极小值得到的拟合键角记为κi，扭转角记为τi，能量函数参数值，记为(c,m,b,d,e,q)；i与步骤(2)中氨基酸的序号一致；N代表氨基酸的总数；

其中，步骤(3).2中确定了参数值的能量函数描述蛋白质的自由能，得到的拟合键角与扭转角描述构建的能量极小构象，简记为拟合结构；最优结构通过将此组新的拟合键角与扭转角构建的拟合结构与步骤(1)的蛋白质实验结构之间的中心碳原子Cα的均方根偏差(RMSD)比较得到；

步骤(4)：模拟细胞环境变化导致固有无序蛋白的构型变化过程；

其中，模拟固有无序蛋白的构型变化过程可采用模拟退火算法，具体为：

细胞环境变化导致固有无序蛋白的构型变化过程是一个非平衡动力学过程，由步骤(3).2输出的拟合结构中的拟合键角与扭转角作为蛋白质的初始状态解，能量函数，即公式(1)，用来计算每一Monte Carlo步的蛋白质能量，对应此Monte Carlo步的新解的接受准则采用Metroplis准则；

Monte Carlo步数根据蛋白质的长度进行合理设置，模拟退火算法中的低温和高温值由实际模拟情况设定；固有无序蛋白从部分结构化到完全无结构化对应着由低温到高温的加热过程，从完全无结构化到部分结构化对应着由高温到低温的冷却过程；

低温，记为T1；高温，记为T2；

步骤(5)：基于步骤(3)的拟合结构，重复步骤(4)对蛋白质进行千次以上的加热-冷却模拟，称为加热-冷却模拟循环，确定具有多态结构的低温值和高温值；

其中，加热-冷却模拟指的是模拟蛋白质从部分结构化到完全无结构化后回到部分结构化的非平衡动力学过程；

步骤(5)中的每一次加热-冷却模拟循环过程为：

步骤(5).1将蛋白质的拟合结构在低温做N0步热平衡；

步骤(5).2用N1步将蛋白质体系线性地从低温加热到高温；

步骤(5).3在高温处对蛋白质进行N2步的平衡热浴，确保蛋白质充分地打开和热平衡；

步骤(5).4接下来的N3步将蛋白质体系线性地从高温冷却回低温；

步骤(5).5最后，在低温做同样的N0步热平衡，使蛋白质冷却充分；

其中，步骤(5)确定具有多态结构的低温和高温值，具体为：

计算加热-冷却模拟循环中的中心碳原子Cα原子RMSD、回旋半径(Rg)和体系能量(E)，观察RMSD和Rg相对于平均值的波动，以低温时RMSD和Rg偏差小于0.1埃为标准确定低温的T1值；一次加热-冷却模拟循环后，初末态构象之间的Cα原子RMSD小于1.0埃，认为蛋白质重新折叠回去；选择不同高温值进行反复测试，统计所有加热-冷却循环后蛋白质折叠回去所占的百分比，比率约为50％的高温值为选定的T2值；

步骤(6)：提纯步骤(5)加热-冷却模拟循环后的蛋白质结构，获得提纯的末态构型数据组，具体为：

经过步骤(5)千次以上加热-冷却模拟循环后，通过体积排斥法提纯加热-冷却循环后的蛋白质结构，即获得提纯的末态构型数据组，具体为：

将步骤(5)所有加热-冷却模拟循环后的蛋白质末态构型输出，以pdb文件的标准格式存储，一个文件存放一个末态构型；采用侧链安装工具PULCHRA软件(http://cssb.biology.gatech.edu/skolnick/files/PULCHRA)或SCWRL4软件(http://dunbrack.fccc.edu/scwrl4/)给所有末态构型安装侧链，得到全原子构型；计算蛋白质分子中任意一对非成键原子的间距，任意原子对间距都大于等于设定值时，视该构型为合理构型而保留；如果存在原子对间距小于设定值，则去掉该构型；

其中，设定值范围在1.2-1.8埃；

步骤(7)：由步骤(6)提纯的末态构型数据组，归纳和分析固有无序蛋白的结构聚类，具体为：

从末态构型数据组中读取Cα原子坐标，计算RMSD、Rg、末端距(Dee)和E；末态能量E由步骤(3)拟合的能量函数，将末态构型的键角与扭转角值直接代入计算；Dee是蛋白质氮末端与碳末端的Cα原子距离；获得所有末态构型的RMSD、Rg、Dee和E后，画柱状分布图或者Rg-Dee、RMSD-Rg、Rg-E二维图或者RMSD-Rg-E、Dee-Rg-E三维图，由这些图展现固有无序蛋白的结构分布的局域性，揭示其结构聚类；

步骤(8)：可视化固有无序蛋白结构聚类的特征构型，具体为：

利用步骤(8)末态构型RMSD、Rg、Dee、E后分布图，找到末态构型的分布区域；其特点是，在分布图中，末态构型的对应点聚集在M个狭小的区间，每个聚集区之间有一定的间隔且分布点相对较多；在每个聚集区随机挑选K个点，找到其对应的结构；用VMD和Jmol为主的蛋白质可视化软件，显示出这些构型的三维结构并进行结构特性分析；

其中，M值为大于1的整数；K值大于等于1的整数；

至此，从步骤(1)到步骤(8)，完成了一种统计计算固有无序蛋白多态结构的方法。

有益效果

一种统计计算固有无序蛋白多态结构的方法，与现有技术及方法相比，具有如下有益效果：

1.直接采用自构建的能量函数进行非平衡动力学模拟，不需要精确描述原子或基团之间相互作用的力场和迭代求解含时牛顿动力学方程，极大地减少计算量；

2.蛋白质折叠和动力学过程计算效率高、速度快，即使对大分子量的蛋白质进行模拟，在普通台式机上几分钟即可完成；

3.非常适合对蛋白质非平衡动力学过程进行大量循环模拟，分析蛋白质折叠和结构的统计学特性；

4.通过大量循环模拟和统计分析，可以揭示固有无序蛋白可能存在的多态构型，解决了实验方法和传统理论方法难以同时捕捉这些多态结构的难题；

5.便于研究固有无序蛋白的构型多样性、潜在的低聚物等，揭示这些蛋白质的细胞毒性和相关疾病的致病机制。

附图说明

图1为一种统计计算固有无序蛋白多态结构的方法流程图；

图2为一种统计计算固有无序蛋白多态结构的方法在具体实施时针对淀粉样-β42肽的流程示意图；

图3为肽淀粉样-β42肽residues 8-40区域的所有提纯后末态构型Dee-Rg分布图；

图4是图3中标示的13个结构聚类对应的特征构型用Jmol可视化的卡通模式图。

具体实施方式

下面结合附图和实施例对本发明的方法作进一步说明。

实施例1

本实施例详细阐述了本发明“一种统计计算固有无序蛋白多态结构的方法”在具体实施时针对一种典型的固有无序蛋白淀粉样-β42肽的计算流程；

图1为一种统计计算固有无序蛋白多态结构的方法的流程图，从图中可以看出，本方法计算过程为：步骤(1)：获取蛋白质实验结构；步骤(2)：采用步骤(1)的蛋白质实验结构中心碳原子坐标数据，计算中心碳原子链接形成的键角和扭转角；步骤(3)：构建蛋白质拟合模型；步骤(4)：模拟细胞环境变化导致固有无序蛋白的构型变化过程；步骤(5)：基于步骤(3)的拟合结构，重复步骤(4)对蛋白质进行千次以上的加热-冷却模拟；步骤(6)：提纯步骤(5)加热-冷却循环后的蛋白质结构；步骤(7)：由步骤(6)提纯的末态构型数据组，归纳和分析固有无序蛋白的结构聚类；步骤(8)：可视化固有无序蛋白结构聚类的特征构型；

图2为本实施例的流程图，从图中可以看出，统计计算固有无序蛋白淀粉样-β42肽多态结构包含如下步骤：

步骤一、获取固有无序蛋白淀粉样-β42肽的实验结构；

一种优选方案是从蛋白质数据银行(PDB)下载淀粉样-β42肽的实验结构1iyt.pdb；

其中，1iyt.pdb是NMR解析的结构，含有10个model的结构信息，本实施例的优选方案是采用model 1的residues 8-40区域统计计算多态结构；

步骤二、读取步骤一实验结构1iyt.pdb中model 1的residues 8-40中心碳原子Cα坐标，计算中心碳原子链接形成的键角和扭转角(κ⁰i,τ⁰i)；

其中，i表示淀粉样-β42肽氮末端到碳末端编号，对应的取值i＝8,9,…,40；

步骤三、构建淀粉样-β42肽residues 8-40区域的拟合模型，具体为：

步骤三.1：观察步骤二输出键角和扭转角的值随编号i的变化，键角和扭转角具有6个常数区域、5个变化区域，从而将residues 8-40区域中心碳原子链相应的分为8-12、13-20、21-27、28-31、32-40五个片段；根据发明内容所述方法步骤(3).1，对每个片段搜索能量函数的极小值，得到每个片段的拟合键角与扭转角及相应能量函数的参数值；

步骤三.2：将步骤三.1中的5个片段顺次拼接在一起，再次搜索能量函数的全局极小，获得若干拟合结构及最优结构；

其中，步骤三.1和步骤三.2中的能量函数为发明内容所述公式(1)；

其中，步骤三.2中确定了参数值的能量函数描述淀粉样-β42肽residues 8-40区域的自由能，得到的拟合键角与扭转角描述构建的能量极小构象，即为淀粉样-β42肽residues 8-40区域的拟合结构；最优结构与步骤一的实验结构的RMSD为1.06埃；

步骤四：模拟细胞环境变化导致淀粉样-β42肽residues 8-40区域的构型变化过程，具体为：

由步骤三.2输出的拟合结构中的拟合键角与扭转角作为淀粉样-β42肽residues 8-40区域的初始状态解，输出的能量函数用来计算每一Monte Carlo步的能量；Monte Carlo步数设置为8×10⁷；模拟退火计算过程中低温取T1＝1×10^-16、1×10^-15、1×10^-14三组测试值，高温取T2＝1×10^-7、1×10^-6、1×10^-5、1×10^-4四组测试值；

步骤五：在每组低温和高温值，基于步骤三的拟合结构，重复步骤四对淀粉样-β42肽residues 8-40区域进行3000-3500次的加热-冷却模拟，输出具有多态结构的加热-冷却模拟后所有末态构型；

在步骤五中每一次加热-冷却模拟循环过程为：

步骤五.1将淀粉样-β42肽residues 8-40的拟合结构在低温做1×10⁷步热平衡；

步骤五.2用2×10⁷步将淀粉样-β42肽residues 8-40线性地从低温加热到高温；

步骤五.3在高温进行2×10⁷步步的平衡热浴，确保淀粉样-β42肽residues 8-40充分地打开和热平衡；

步骤五.4接下来的2×10⁷步将淀粉样-β42肽residues 8-40线性地从高温冷却回低温；

步骤五.5最后，在低温做同样的1×10⁷步热平衡，使淀粉样-β42肽residues 8-40冷却充分；

其中，在步骤五需要确定具有多态结构的低温和高温值，具体为：

计算加热-冷却模拟循环中的淀粉样-β42肽residues 8-40中心碳原子Cα原子RMSD和回旋半径Rg，观察RMSD和Rg相对于平均值的波动，发现低温T1＝1×10^-15时RMSD和Rg偏差小于0.1埃；统计所有加热-冷却循环后蛋白质折叠回去所占的百分比，发现高温T2＝1×10^-6时折叠回去百分比约为50％；从而，确定具有多态结构的低温和高温为：T1＝1×10^-15、T2＝1×10^-6；

步骤六：提纯步骤五淀粉样-β42肽residues 8-40加热-冷却循环后后具有多态结构的末态构型，具体为：

将步骤五在低温T1＝1×10^-15和高温T2＝1×10^-6时加热-冷却模拟后的所有末态构型，采用SCWRL4安装侧链得到全原子构型；采用1.6埃的截断值进行原子体积排斥计算，剔除不合理构型，得到提纯后的3000个末态构型；

步骤七：由步骤六提纯的3000个末态构型，采用Dee-Rg分布图给淀粉样-β42肽residues 8-40的末态构型分类，具体为：

根据发明内容所述方法步骤(7)，从3000个末态构型中读取Cα原子坐标，计算RMSD、Rg、末端距Dee；图3是3000个提纯后末态构型的Dee-Rg分布图，图中横坐标是回旋半径Rg，纵坐标是末端距Dee；图中用椭圆框圈出了的主要聚集区域，椭圆框顺序以数字1-13标记，它们对应着13个主要聚集区域；每个区域的末态构型划分为一个结构聚类，1-13聚集区域对应的末端距Dee值和回旋半径Rg值的范围分布为：聚集区1为Dee∈[37,38]、Rg∈[11.75,11.909]，聚集区2为Dee∈[35.2,37]、Rg∈[11.2,11.6]，聚集区3为Dee∈[32,34]、Rg∈[11.155,11.4]，聚集区4为Dee∈[28,29.4]、Rg∈[10.4,10.6]，聚集区5为Dee∈[25,26.5]、Rg∈[9.5,9.61]，聚集区6为Dee∈[22.5,23.2]、Rg∈[10.4,10.6]，聚集区7为Dee∈[20.6,22.2]、Rg∈[9.6,9.7]，聚集区8为Dee∈[19,22.5]、Rg∈[9.3,9.52]，聚集区9为Dee∈[19.6,20.9]、Rg∈[9.9,10.1]，聚集区10为Dee∈[15,18]、Rg∈[8.92,9.3]，聚集区11为Dee∈[20,22.4]、Rg∈[8.15,8.4]，聚集区12为Dee∈[16,17]、Rg∈[8.05,8.16]，聚集区13为Dee∈[10.9,11.8]、Rg∈[7.85,7.95]；

步骤八：采用Jmol软件可视化淀粉样-β42肽residues 8-40结构聚类的特征构型，具体为：

在步骤七末态构型Dee-Rg分布图，确定了淀粉样-β42肽residues 8-40的13个结构聚类的分布区域；可以发现，每个结构聚类区域包含许多末态构型的对应点，每个聚集区之间有一定的间隔且分布点相对较多；每个结构聚类中对应的末态构型Dee、Rg相差很小，这表明它们的结构是相似的；不同结构聚类之间具有较大的构型差异；图4是在结构聚类中随机挑选1个点，找到其对应的pdb结构，然后采用Jmol软件画出该结构卡通模式图；图4中结构图编号1-13与图3一致，比较不同的结构聚类的特征构型，可以发现它们的螺旋、环区域存在的差异。

以上所述为本发明的较佳实施例而已，本发明不应该局限于该实施例和附图所公开的内容。凡是不脱离本发明所公开的精神下完成的等效或修改，都落入本发明保护的范围。