【2.5.1.5】并行计算蛋白蛋白结合自由能

一、建立起始结构并运行模拟以获得平衡的系统。

我们将在此模拟中建模的系统是人类H-Ras蛋白与C-Raf1的Ras结合域(Ras-Raf)之间的复合体,它是信号转导级联的核心。 这是RAS-RAF复合系统的部分平衡的,预先准备的pdb文件。

ras-raf.pdb

该结构包含ras和raf蛋白,以及生理必需的GTP核苷酸,如下图所示:

有关如何构建起始结构和运行仿真以获得平衡系统的指导,请参阅第1节第2节

使用MMPBSA.py计算绑定自由能的重要文件是拓扑文件和mdcrd文件( ras-raf_top_mdcrd.tgz

二、并行计算Ras-Raf的结合自由能。

现在,我们将计算复合物,受体和配体的相互作用能和溶剂化自由能,并对结果求平均值,以获得结合自由能的估计值。请注意,在本教程的这一部分中,我们将不会计算熵对绑定的贡献,因此严格来说,我们的结果将不是真正的自由能,但可用于与类似系统进行比较。有关使用正态分析(Nmode)来计算系统的熵贡献或取消注释常规部分中最后一行以使用AMBER中的ptraj模块执行准谐波熵计算的示例,请参见第5节

我们将同时使用MM-GBSA和MM-PBSA方法进行结合能计算,以与从3.1节 获得的结果进行比较。 MMPBSA.py.MPI通过为每个线程(进程)分配相等数量的帧来并行化计算。因此,当处理的帧数是启动的线程数的倍数时,它最有效地运行。但是,这不是必需的。如果帧数不是线程数的倍数,则“剩余”帧将在启动的线程数的子集中平均分配。例如,在3个线程上运行50帧将导致2个线程计算17个帧,而最后一个线程仅计算16个帧。因此,第三个线程将必须等待前两个线程完成计算才能继续。因此,例如5个线程将​​是一个明智的选择(因为每个线程占用10帧)。但是,线程数不能超过已处理的帧数,否则MMPBSA.py.MPI将终止并显示一条错误消息。 MMPBSA.py.MPI的输入文件与MMPBSA.py的输入文件相同:

mmpbsa.in

Input file for running PB and GB
&general
   endframe=50, verbose=1,
#  entropy=1,
/
&gb
  igb=2, saltcon=0.100
/
&pb
  istrng=0.100,
/

MMPBSA.py.MPI的输入文件被设计为类似于AMBER的sander模块中使用的mdin文件的设置。 每个名称列表的开头均以&符号表示,后跟名称列表的名称。 此外,反斜杠(/)或’&end’可用于结束名称列表。 有关所有变量的完整列表,请参见amber手册。 该输入文件分为三个名称列表:general,pb和gb。 通用名称列表旨在指定并非特定于计算的特定部分而是特定于所有部分的变量。 在此设置中,我们已将RAS定义为受体,将RAF定义为配体。 ‘endframe’变量设置mdcrd的停止帧。 “&gb”和“&pb”名称列表标记使脚本知道使用这些名称列表中定义的给定值执行MM-GBSA和MM-PBSA计算。

mpirun -np 4 $AMBERHOME/bin/MMPBSA.py.MPI -O -i mmpbsa.in -o FINAL_RESULTS_MMPBSA.dat -sp ras-raf_solvated.prmtop -cp ras-raf.prmtop -rp ras.prmtop -lp raf.prmtop -y *.mdcrd > progress.log 2>&1

或通过将脚本提交到排队系统(例如PBS)

qsub parallel.job

parallel.job

#!/bin/sh
#PBS -N rasraf_parallel
#PBS -o parallel.out
#PBS -e parallel.err
#PBS -m abe
#PBS -M email@address.com
#PBS -q brute
#PBS -l nodes=1:node:ppn=4
#PBS -l pmem=900mb

cd $PBS_O_WORKDIR

mpirun -np 4 MMPBSA.py.MPI -O -i mmpbsa.in -o FINAL_RESULTS_MMPBSA.dat -sp ras-raf_solvated.prmtop -cp ras-raf.prmtop -rp ras.prmtop -lp raf.prmtop -y bigprod.mdcrd > progress.log 2>&1

这会将计算进度打印到文件 progress.log 中。 计算过程中的所有错误也将被打印到该文件中(这是2>&1的目的)。 最后,一旦计算完成,将打印计时,显示脚本每个步骤所花费的时间。

progress.log

MMPBSA.py.MPI being run on 4 processors
ptraj found! Using /scr/arwen_3/swails/i686/amber11/exe/ptraj
sander found! Using /scr/arwen_3/swails/i686/amber11/exe/sander

Preparing trajectories with ptraj...
50 frames were read in and processed by ptraj for use in calculation.

Starting calculations

Starting gb calculation...

Starting pb calculation...


Calculations complete. Writing output file(s)...
Timing:
Processing Trajectories With Ptraj:       0.126 min.
Total GB Calculation Time (sander):       4.782 min.
Total PB Calculation Time (sander):      28.407 min.
Output File Writing Time:                 0.053 min.

Total Time Taken:                        33.379 min.


MMPBSA Finished. Thank you for using. Please send any bugs/suggestions/comments to mmpbsa.amber@gmail.com

命令行上的’-y * .mdcrd’告诉脚本读取工作目录中所有以’.mdcrd’结尾的文件,并将它们用作要分析的轨迹。

将’keep_files’设置为其默认值1,以下是所有输出文件:Parallel_output.tgz。

该脚本使用ptraj创建三个未溶剂化的mdcrd文件(复合物,受体和配体),这是在GB和PB计算期间分析的坐标。 * .mdout文件包含所有指定帧的能量。 如果使用ptraj执行准谐波熵计算,则会将平均pdb文件创建为平均结构。 除最终输出文件FINAL_RESULTS_MMPBSA.dat以外,由MMPBSA.py.MPI创建的所有文件都应以前缀’_MMPBSA_开头。

FINAL_RESULTS_MMPBSA.dat

| Run on Sun Feb 14 19:10:43 EST 2010

|Input file:
|--------------------------------------------------------------
|Input file for running PB and GB in serial
|&general
|   endframe=50, verbose=1,
|   mpi_cmd='mpirun -np 3', nproc=3
|/
|&gb
|  igb=2, saltcon=0.100
|/
|&pb
|  istrng=0.100, 
|/
|--------------------------------------------------------------
|Solvated complex topology file:  ras-raf_solvated.prmtop
|Complex topology file:           ras-raf.prmtop
|Receptor topology file:          ras.prmtop
|Ligand topology file:            raf.prmtop
|Initial mdcrd(s):                bigprod.mdcrd
|
|Best guess for receptor mask:   ":1-166"
|Best guess for  ligand  mask:   ":167-242"

|Calculations performed using 50 frames.
|Poisson Boltzmann calculations performed using internal PBSA solver in sander.
|
|All units are reported in kcal/mole.
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------

GENERALIZED BORN:

Complex:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -1863.7944               17.1704              2.4283
EEL                     -17200.7297               75.9366             10.7391
EGB                      -3249.6511               65.2075              9.2217
ESURF                       91.3565                1.3938              0.1971

G gas                   -19064.5240               77.8536             11.0102
G solv                   -3158.2946               65.2224              9.2238

TOTAL                   -22222.8186               51.0216              7.2155


Receptor:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -1268.1888               14.2342              2.0130
EEL                     -11557.0773               71.7127             10.1417
EGB                      -2532.0669               57.7003              8.1600
ESURF                       64.2843                1.1143              0.1576


G gas                   -12825.2661               73.1118             10.3396
G solv                   -2467.7826               57.7110              8.1616

TOTAL                   -15293.0487               35.3527              4.9996


Ligand:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                   -529.3090                9.4198              1.3322
EEL                      -4684.4720               36.1449              5.1117
EGB                      -1688.9631               26.5353              3.7527
ESURF                       37.0493                0.6185              0.0875


G gas                    -5213.7811               37.3522              5.2824
G solv                   -1651.9138               26.5425              3.7537

TOTAL                    -6865.6949               25.8878              3.6611


Differences (Complex - Receptor - Ligand):
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                    -66.2966                4.2751              0.6046
EEL                       -959.1803               34.9190              4.9383
EGB                        971.3789               33.0497              4.6739
ESURF                       -9.9770                0.3759              0.0532


DELTA G gas              -1025.4769               35.1797              4.9752
DELTA G solv               961.4018               33.0518              4.6742


 DELTA G binding =        -64.0750     +/-      6.3729                 0.9013
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------

POISSON BOLTZMANN:

Complex:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -1863.7944               17.1704              2.4283
EEL                     -17200.7297               75.9366             10.7391
EPB                      -3207.7160               66.4023              9.3907
ECAVITY                     67.8762                0.7818              0.1106

G gas                   -19064.5240             6061.1875            857.1813
G solv                   -3139.8399               66.4069              9.3914

TOTAL                    -7686.8660               52.5400              7.4303


Receptor:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                  -1268.1888               14.2342              2.0130
EEL                     -11557.0773               71.7127             10.1417
EPB                      -2483.7242               56.4551              7.9840
ECAVITY                     47.1495                0.4737              0.0670

G gas                   -12825.2661             5345.3320            755.9441
G solv                   -2436.5747               56.4571              7.9842

TOTAL                    -5250.2060               38.5188              5.4474


Ligand:
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                   -529.3090                9.4198              1.3322
EEL                      -4684.4720               36.1449              5.1117
EPB                      -1670.4169               27.6694              3.9131
ECAVITY                     28.0328                0.4133              0.0584

G gas                    -5213.7811             1395.1865            197.3092
G solv                   -1642.3841               27.6725              3.9135

TOTAL                    -2350.3020               25.1197              3.5525


Differences (Complex - Receptor - Ligand):
Energy Component            Average              Std. Dev.   Std. Err. of Mean
-------------------------------------------------------------------------------
VDWAALS                    -66.2966                4.2751              0.6046
EEL                       -959.1803               34.9190              4.9383
EPB                        946.4251               34.5128              4.8808
ECAVITY                     -7.3062                0.3004              0.0425

DELTA G gas              -1025.4769             1237.6138            175.0250
DELTA G solv               939.1189               34.5141              4.8810



 DELTA G binding =        -86.3579     +/-      8.3264                 1.1775
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------

WARNINGS:
igb=2 should be used with mbondi2 pbradii set. Yours are modified Bondi radii (mbondi)

统计文件的开头包括日期/时间,输入文件的副本,配体和受体掩膜的最佳猜测,脚本使用的文件,分析的帧数(frames)以及哪个PB solver(如果有)被使用了。统计信息文件的其余部分包括GB的所有平均能量,标准偏差和均值的标准误差,其后是PB。在每个部分之后,结合的ΔG随误差值一起给出。该文件中不同术语的含义如下:

VDWAALS = MM的van der Waals贡献。
EEL =由MM力场计算出的静电能。
EPB / EGB =分别由PB或GB计算得出的对溶剂化自由能的静电贡献。
ECAVITY =通过经验模型计算得出的对溶剂化自由能的非极性贡献。
DELTA G binding=根据上述条件计算得出的最终估计结合自由能。 (千卡/摩尔)

请注意,尚未报告总气相能,因为使用单轨方法,受体和配体的键合项的值应恰好抵消了络合物的键合项的值。如果能量没有抵消到允许的公差(0.001 kcal / mol)之内,则会出现错误消息。

人们通常期望找到一种非常有利的静电能和一种不利的溶剂化自由能。这象征着人们必须使用的能量来使结合颗粒去溶剂化并对齐其结合界面。

从负的总结合自由能为-86.36 kcal / mol,我们可以清楚地看到,这是纯水中的一种有利的蛋白质-蛋白质复合物,但请记住,由于我们没有估算(不利于)熵对结合的贡献。请注意,GB方法给出的结合能略低,但仍表明这是一种有利的结合态。

参考资料

个人公众号,比较懒,很少更新,可以在上面提问题,如果回复不及时,可发邮件给我: tiehan@sina.cn

Sam avatar
About Sam
专注生物信息 专注转化医学