一种分析蛋白质侧链构象含时动力学演化的方法与流程

本发明涉及分析蛋白质折叠或其他动力学过程中空间结构的变化，尤其涉及一种分析蛋白质侧链构象含时动力学演化的方法，属于蛋白质结构与动力学研究领域。

背景技术：

蛋白质是最重要的生物大分子之一，它在生物体内承担着物质输运、信号传递和能量供给等生物学功能。蛋白质的生物学功能是由其自然折叠而成的空间结构决定的。实验上，采用x射线衍射、多维核磁共振等技术可以测定蛋白质的结构。理论上，人们发展了分子动力学技术研究蛋白质折叠和动力学过程。它通过计算蛋白质分子中原子的运动轨迹、原子间的相互作用等，在原子尺度上阐明蛋白质实现生物学功能过程中空间结构演化的微观动力学机制。然而，蛋白质分子的组成与结构非常复杂，其折叠和动力学过程涉及大量原子的集体运动，发展一些结构分析方法非常必要。目前已发展了一些理论方法，以及空间结构可视化分析软件，如cn3d、vmd、jmol等。这些方法和工具在展示蛋白质分子三维空间构象、揭示结构特征、分析原子或基团间相互作用等方面非常有用。但是，它们是对某一时刻或某一帧的空间结构进行分析，缺乏展现蛋白质分子空间结构含时演化的能力。

由于组成蛋白质的各种氨基酸主要区别是侧链基团，每个蛋白质所具有的独特空间结构和功能与侧链紧构型密相关。在蛋白质实现其生物学功能的过程中，相关侧链构型的伴随变化(如旋转、扭转、弯曲等)是一个重要的方面。因而，在蛋白质折叠或其他动力学过程模拟中，分析蛋白质的侧链构象含时演化显得非常重要，目前这方面的研究还很少。蛋白质的空间结构演化是主链和侧链的协同运动，目前的结构分析工具没有将主链和侧链运动分离，难以观测特定侧链的构象演化，还需要更有效的分析手段和方法。

技术实现要素：

本发明的目的是针对目前还没有分析蛋白质折叠和其他动力学过程中侧链结构变化的技术现状，提出了一种分析蛋白质侧链构象含时动力学演化的方法。

本发明所提方法通过计算蛋白质运动过程中氨基酸残基侧链在每一个时刻的相对扭转角，获得氨基酸残基侧链扭转为主的重要信息。

为实现上述目的，一种分析蛋白质侧链构象含时动力学演化的方法，步骤如下：

步骤1：获取蛋白质初始结构；

其中，蛋白质初始结构可以从蛋白质数据银行(pdb，http://www.rcsb.org)获取，也可以通过蛋白质结构建模的方法与工具获取；

步骤2：采用步骤1的蛋白质初始结构进行含时分子动力学模拟，得到蛋白质运动轨迹文件，具体为：

以步骤1的蛋白质初始结构，利用分子动力学模拟软件，选择合适的力场，对蛋白质进行含时分子动力学模拟，得到蛋白质中各原子的运动轨迹信息，保存运动轨迹文件；

其中，力场的一种优选方案是全原子力场和联合原子力场；

其中，模拟积分步长取δt飞秒，对整个系统进行t纳秒的分子动力学模拟；每隔w个积分步输出一次蛋白质的结构，即每δt×w飞秒记录一次蛋白质结构的信息；运动轨迹文件共保存k＝(t×10⁶)/(δt×w)个时刻的结构，依次编号为n＝1,2,3,…,k；

其中，δt、w、k值为大于1的正整数；

步骤3：从步骤2的蛋白质运动轨迹文件，提取氨基酸残基侧链，计算第n时刻的氨基酸残基侧链的相对扭转角，具体为：

从步骤2得到的蛋白质运动轨迹文件，获得第n时刻蛋白质的结构，输出坐标数据文件；从坐标数据文件中，提取两个待考察氨基酸残基中主链碳原子、侧链第一个碳原子的坐标；计算这四个碳原子形成的二面角，即为两个待考察氨基酸残基侧链的相对扭转角；

其中，第n时刻对应步骤2中编号n等于1的时刻；

其中，主链碳原子、侧链第一个碳原子记为cα、cβ，它们从蛋白质的氮末端到碳末端根据氨基酸残基顺序编号；两个待考察氨基酸残基的编号记为i和j，i和j为1到m的正整数，m是蛋白质中氨基酸残基总数；i表示编号靠前的氨基酸残基，j表示编号靠后的氨基酸残基，j大于i；j与i的差等于1表示两个氨基酸残基相邻，j与i的差大于1表示两个氨基酸残基非相邻；

其中，第i个氨基酸残基的主链碳原子、侧链第一个碳原子记为cⁱα、cⁱβ；它们坐标分别记为rⁱα、rⁱβ；第j个氨基酸残基的主链碳原子、侧链第一个碳原子记为c^jα、c^jβ；它们坐标分别记为r^jα、r^jβ；

两个待考察氨基酸残基侧链的相对扭转角，具体通过如下步骤计算：

步骤3.1计算第i个氨基酸残基的cⁱα原子指向第j个氨基酸残基的c^jα原子的单位矢量x^ijα，表述为如下公式(1)：

步骤3.2计算第i个氨基酸残基的cⁱα原子指向侧链cⁱβ原子的单位矢量aⁱβ，表述为如下公式(2)：

步骤3.3计算步骤3.2的单位矢量aⁱβ在步骤3.1的单位矢量x^ijα法平面上的单位投影矢量aⁱ，表述为如下公式(3)：

步骤3.4计算第j个氨基酸残基的c^jα原子指向侧链c^jβ原子的单位矢量b^jβ，表述为如下公式(4)：

步骤3.5计算步骤3.4的单位矢量b^jβ在步骤3.1的单位矢量x^ijα法平面上的单位投影矢量b^j，表述为如下公式(5)：

步骤3.6由步骤3.3和步骤3.5的单位投影矢量aⁱ、b^j，计算cⁱβ-cⁱα-c^jα-c^jβ原子形成的二面角φ，即第i和j个氨基酸残基侧链的相对扭转角，表述为如下公式(6)：

cosφ＝aⁱ·b^j(6)

步骤4：计算第一时刻后所有时刻的氨基酸残基侧链的相对扭转角，给出相对扭转角随时间的演化，即：重复步骤3，并且n对应n＝2,3,…,k，具体为：

从步骤2得到的蛋白质运动轨迹文件，获得第一时刻后所有时刻的蛋白质结构，输出相应的坐标数据文件；采用与步骤3相同的方法，计算两个待考察氨基酸残基侧链在这些时刻的相对扭转角，输出保存这些时刻相对扭转角的值；以时间为横轴，相对扭转角的值为纵轴，画出从第一时刻到最终时刻氨基酸残基侧链的相对扭转角的变化图；

其中，两个待考察氨基酸残基与步骤3中的提取过程一致；

其中，第一时刻后所有时刻对应步骤2描述中的n＝2,3,…,k时刻；

至此，从步骤1到步骤4，完成了一种分析蛋白质侧链构象含时动力学演化的方法。

有益效果

一种分析蛋白质侧链构象含时动力学演化的方法，与现有技术及方法相比，具有如下有益效果：

(1)能够分析蛋白质折叠和其他动力学过程中蛋白质结构的含时演化；

(2)能够将蛋白质含时动力学过程中主链和侧链运动分离，分析特定侧链的构象演化；

(3)不仅能够分析非相邻的氨基酸残基侧链的含时动力学演化，也可分析所有相邻的氨基酸残基侧链的含时动力学演化；

(4)对研究蛋白质折叠过程中二级结构的形成、蛋白质生物学功能实现与构型变化内在联系等具有基础和应用价值；

(5)便于在与蛋白质折叠和动力学、蛋白质结构和功能相关的生命科学及生物医药学领域广泛应用。

附图说明

图1为一种分析蛋白质侧链构象含时动力学演化的方法流程图；

图2为原癌基因c-myc蛋白质的含时动力学演化过程中leu951与ala955侧链相对扭转角变化计算的流程示意图；

图3为原癌基因c-myc蛋白质中leu951与ala955的侧链相对扭转角随时间的演化图；

图4是原癌基因c-myc蛋白质中所有相邻氨基酸残基侧链的相对扭转角在n＝[300-1300]时刻的演化图。

具体实施方式

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

实施例1

本实施例详细阐述了本发明“一种分析蛋白质侧链构象含时动力学演化的方法”在具体实施时针对原癌基因c-myc蛋白质的含时动力学演化过程中leu951与ala955侧链相对扭转角变化计算的流程。

图1是一种分析蛋白质侧链构象含时动力学演化的方法的流程图。从图中可以看出，本方法包含的过程为：步骤(1)：获取蛋白质初始结构；步骤(2)：采用步骤(1)的蛋白质初始结构进行含时分子动力学模拟，保存蛋白质运动轨迹文件；步骤(3)：计算第一个时刻的氨基酸残基侧链的相对扭转角；步骤(4)：计算其他时刻的氨基酸残基侧链的相对扭转角，给出相对扭转角随时间的演化；

图2为本实施例的流程图，从图中可以看出，原癌基因c-myc蛋白质的含时动力学演化过程中leu951与ala955侧链相对扭转角变化计算包含如下步骤：

步骤一、从蛋白质数据银行(pdb，http://www.rcsb.org)下载代码为1nkp的结构数据文件，提取原癌基因c-myc蛋白质的结构，保存为1nkp.pdb；

步骤二、以步骤一的1nkp.pdb为初始结构，利用gromacs4.6.3软件，采用gromos53a6力场，对原癌基因c-myc蛋白质进行含时分子动力学模拟，得到蛋白质中各原子的运动轨迹，保存为1nkp_md.trr；

其中，蛋白质置于长方体水盒子中，蛋白质到水盒子壁的距离设为2.0纳米；水分子采用spc模型；nacl浓度设为0.15摩尔/升；模拟中，采用周期性边界条件；温度设为290开，压强取1个标准大气压，平衡时采用berendsen方法控温控压，分子动力学计算时采用vrescale方法控温、parrinelo-rahman方法控压；长程相互作用的截断半径取0.9埃；积分步长取2飞秒，对整个系统进行50纳秒的分子动力学模拟；

其中，每20皮秒记录一次蛋白质结构的信息，1nkp_md.trr文件中共保存2500个时刻的结构，依次编号为n＝1,2,3,…,2500，对应着时间t＝20,40,60,…,50000皮秒时蛋白质的结构；

步骤三、计算第1个时刻leu951与ala955侧链的相对扭转角，具体为：

从步骤二得到的蛋白质中各原子的运动轨迹1nkp_md.trr文件，获得n＝1时刻蛋白质的结构，提取非相邻的两个氨基酸残基leu951与ala955中主链cα、侧链cβ原子的坐标r⁹⁵¹α、r⁹⁵¹β、r⁹⁵⁵α、r⁹⁵⁵β；按照发明内容步骤(3)中步骤(3).1到步骤(3).5所述，采用公式(1)到公式(6)，计算出n＝1时刻两个氨基酸残基leu951与ala955侧链的相对扭转角，记为

其中，的下标与步骤二中n的编号一致；

步骤四、计算第2到2500时刻leu951与ala955侧链的相对扭转角，给出相对扭转角随时间的演化，具体为：

从步骤二得到的蛋白质中各原子的运动轨迹1nkp_md.trr文件，获得n＝2,3,…,2500时刻蛋白质的结构，提取这些时刻氨基酸残基leu951与ala955中主链cα、侧链cβ原子的坐标；在每一个时刻，采用发明内容步骤(3)中步骤(3).1到步骤(3).5所述的公式(1)到公式(6)，计算出两个氨基酸残基leu951与ala955侧链的相对扭转角，其值按照时刻编号顺序记为以时刻编号n为横轴，相对扭转角的值为纵轴，画出氨基酸残基leu951与ala955侧链的相对扭转角的变化图；

图3是原癌基因c-myc蛋白质中leu951与ala955的侧链相对扭转角随时间的演化图；图中横坐标xaxis是时刻编号n，纵坐标yaxis是leu951与ala955的侧链相对扭转角；图中实线是相应时间区间中相对扭转角的平均值，虚线是该时间区间中相对扭转角的平均波动；bend表示弯曲构型，turn表示转角构型；从图中相对扭转角随时间的演化可以分析出：n＝[1,400]区间的相对扭转角平均值约为0.52，n＝[400,500]区间的相对扭转角平均值约为-0.29，n＝[500,700]区间的相对扭转角平均值约为0.46，n＝[700,800]区间的相对扭转角平均值约为0.30和-0.45；这个结果反映出，在该动力学过程中，leu951与ala955的侧链开始处于较小的相对扭转，到n＝[400,500]时会出现小的反向相对扭转，n＝[500,700]会变回与[1,400]相近的扭转角，[700,800]时存在正反扭转之间快速变换；

至此，从步骤一到步骤四，完成了原癌基因c-myc蛋白质的含时动力学演化过程中leu951与ala955侧链相对扭转角变化计算。

实施例2

本实施例按照本发明“一种分析蛋白质侧链构象含时动力学演化的方法”的步骤和实例例1所述流程，阐述原癌基因c-myc蛋白质的含时动力学演化过程中所有相邻氨基酸残基侧链的相对扭转角变化计算及其结果。

原癌基因c-myc蛋白质的含时动力学演化过程中所有相邻氨基酸残基侧链的相对扭转角变化计算，步骤a、b与实施例1步骤一、二相同；在进行氨基酸残基侧链相对扭转角随时间的演化计算时，需依次选取原癌基因c-myc蛋白质的氨基酸残基对，方法是从第1个氨基酸残基开始按顺序选择：1与2、2与3、3与4、…、m-1与m，遇到甘氨酸残基跳过不计，m是蛋白质中氨基酸残基总数；每个氨基酸残基对的侧链相对扭转角计算与实施例1步骤三、四相同；最终给出的是，所有时刻、所有相邻的氨基酸残基侧链相对扭转角；

图4是原癌基因c-myc蛋白质中所有相邻氨基酸残基侧链的相对扭转角在n＝[300-1300]时刻的演化图；图中横坐标xaxis是时刻编号n，纵坐标yaxis是的氨基酸残基编号；一种优选的办法是，用相对扭转角值与颜色的对应表示角度值变化，蓝色到红色依次对应着-π到+π弧度；图中显示出，在n＝[560,800]区间，959-960、960-961、961-962的氨基酸残基对的侧链相对扭转角会明显增大；在n＝[800,900]区间，955-956氨基酸残基对的侧链相对扭转角也会变大；

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