上期四川大学李建宗博士为大家分享了》
。。。。

今天，，，，将继续为大家介绍如何在舜王生物云计算平台上利用Amber完成蛋白-小分子体系的动力学模拟
，，以及利用MMGBSA计算小分子与蛋白质（Abl和伊马替尼）之间的结合自由能
。。。。

Amber是美国加州大学Peter Kollman等开发的一款著名的分子动力学模拟软件包。。
。Amber主要适用于蛋白质
，
，，，小分子和多糖等生物分子体系的模拟。。。
。

建议现在自己电脑上安装一个UCSF Chimera小程序，
，
，，用于处理分子模拟所需的结构文件。。该软件学术免费，
，而且具备和PyMOL类似的图形显示功能。。获取地址：

https://www.cgl.ucsf.edu/chimera/download.html

• 获取晶体结构结构复合物

本教程旨在利用Amber在舜王生物云平台上完成蛋白质Abl和靶向药物伊马替尼体系的动力学模拟，，，并且计算他们之间的结合自由能。。幸运的是PDB数据库中包含了Abl和伊马替尼的复合物的晶体结构，，编号为6E4F。。通过UCSF Chimera的fetch工具可以方便的获取晶体结构。。。。在File—Fetch by ID中的PDB后输入编号6E4F即可下载复合物晶体结构，
，，，并将其保存为6E4F.pdb到自定义工作路径中
。。。

 Abl和伊马替尼的相互作用细节

• 处理蛋白质

 下载复合物晶体结构以后
，，，，先对蛋白质进行处理
，，，，删除复合物中所有的除开蛋白质外的原子
。
。在Residue中选择所有非标准残基，，然后通过Actions—Atoms/bonds—delete进行删除。。
。。将蛋白质保存为protein.pdb文件。。。。

• 处理配体

复合物中伊马替尼的名称为HRA，，，重新打开6E4F.pdb文件，，然后选中HRA，，在Select菜单中用Invert选中其余分子，
，，
，并且删除，，，只保留伊马替尼的结构
。。。如下所示。。。

然后在Tools中找到Structure Editing分别对小分子进行加氢（Add H）和加电荷（Add charge）处理,电荷选择AM1-BCC。
。并将处理后的小分子保存为ligand.mol2文件
。。

登录舜王生物云计算平台
，
，，将处理好的结构文件上传至舜王生物云计算平台，，然后在提交作业中选择图形界面提交或者命令界面提交。
。。。舜王生物云可快速帮你安装诸如Amber在内的模拟软件，，，免去了安装和环境配置等繁琐步骤，，
，
，而且上传输入文件和下载计算结果文件都是图形操作，，，直观且方便
。。
。

Amber支持GPU加速，，在NVIDIA GPU 上的运行速度比仅使用 CPU 的系统快 15 倍
，，，因此，，在提交作业的时候注意选择GPU计算区，，这里的GPU计算资源十分丰富。。。。对于不需要GPU的作业，，可以选择通用算区节省成本
。。

通过图形界面提交任务，，，然后启动计算节点后，，
，，进入到我们刚才上传的文件夹下，，可以看到如下所示的操作系统界面，，舜王生物云计算平台目前使用的是Centos 7系统，
，鼠标右键选择Open in Terminal然后就可以开始进行动力学模拟和结合自由能计算了。。
。

本教程展示在舜王生物云计算平台上进行动力学模拟的具体细节
，
，计算量大的任务请参考本教程第04小结进行任务提交。。。。

利用Antechamber处理配体，，
，，并且利用parmchk2检查转换为Amber识别的参数结构文是否正确，，命令如下：

antechamber -i ligand.mol2 -fi mol2 -o ligand.ante.mol2 -fo mol2 
parmchk2 -i ligand.ante.mol2 -f mol2 -o ligand.frcmod

其中-i指定输入文件，，，，-fo或则-f指定输入文件格式
，
，，，-o指定输出文件，，，，-fo指定输出文件格式。。。检查是否有ligand.ante.mol2和ligand.frcmod文件生成。
。利用Chimera打开ligand.ante.mol2，
，可观测到如下参数化后的结构。。。。

然后利用tleap模块生成配体
、
、、蛋白以及复合物的拓扑文件和模拟初始文件。。首先调用tleap

tleap

然后在tleap窗口中分别加载小分子、、、蛋白质和溶剂所需要的力场参数

source leaprc.protein.ff14SBsource leaprc.gaffsource leaprc.water.tip3p

其中protein.14SB指定蛋白质力场，，gaff指定小分子力场，，
，，tip3p指定水分子模型。。。。

然后生成配体初始文件ligand.prmtop和ligand.inpcrd
，，
，，并且检测配体结构完整性
。。

LIG = loadmol2 ligand.ante.mol2check LIGloadamberparams ligand.frcmodsaveoff LIG ligand.libsaveamberparm LIG ligand.prmtop ligand.inpcrd

生成蛋白初始文件receptor.prmtop和receptor.inpcrd

REC = loadpdb rec.pdbcheck RECsaveamberparm REC receptor.prmtopreceptor.inpcrd

生成蛋白-配体复合物初始文件com_solvated.prmtop和com_solvated.inpcrd

COM = combine {REC, LIG}saveamberparm COM com.prmtop com.inpcrdcharge COMsolvatebox COM TIP3PBOX 12.0saveamberparm COM com_solvated.prmtopcom_solvated.inpcrd

tleap会有类似提示

Building topology.Building atom parameters.Building bond parameters.Building angle parameters.Building proper torsion parameters.Building improper torsion parameters.total943 improper torsions appliedBuilding H-Bond parameters.Incorporating Non-Bonded adjustments.Not Marking per-residue atom chain types.Marking per-residue atom chain types.(Residues lacking connect0/connect1 -these don't have chain types marked:res total affected

完成上述操作后，，没有错误提示，，，，则可退出tleap

quit

如果成功生成下面这些这些文件，，，则可以开始下一步的能量优化了。。

• 能量最小化

动力学模拟起始需要对体系进行初始能量最小化，，以消除不合理原子间的碰撞和不规范的作用力，，，命令如下：

pmemd.cuda -O -i min.in -o min.out -p com_solvated.prmtop -c com_solvated.inpcrd -r min.rst -ref com_solvated.inpcrd

pmemd.cuda,调用pmemd的GPU版本,这里-i指定参数文件min.in, 对体系的平衡进行了限定，，，，参考官网指南具体min.in内容

minimise com&cntrlimin=1,maxcyc=2000,ncyc=1000,cut=8.0,ntb=1,ntc=2,ntf=2,ntpr=100,ntr=1, restraintmask=':1-288',restraint_wt=2.0/

imin=1 指定模拟步骤为能量最小化

maxcyc表明采用共轭梯度算法
，
，，，且指定最大循环数，，，，

ncyc=1000表明采用最陡下降算法

ntpr=100表明保存模拟信息间隔

restraintmask指定限制氨基酸范围

restraint_wt=2.0表示限制力的大小

• 升温模拟

对体系进行升温，
，温度从0K开始到300K结束，
，命令如下
，，并且检测是否有heat.rst和 heat.mdcrd文件生成。。

pmemd.cuda -O -i heat.in -o heat.out -p com_solvated.prmtop -c min.rst -r heat.rst -x heat.mdcrd -ref min.rst

heat.in参数内容

heat ras-raf&cntrlimin=0,irest=0,ntx=1,nstlim=25000,dt=0.002,ntc=2,ntf=2,cut=8.0, ntb=1,ntpr=500, ntwx=500,ntt=3, gamma_ln=2.0,tempi=0.0, temp0=300.0, ig=-1,ntr=1, restraintmask=':1-242',restraint_wt=2.0,nmropt=1/&wt TYPE='TEMP0', istep1=0, istep2=25000,value1=0.1, value2=300.0, /&wt TYPE='END' /

• 密度平衡

对体系进行密度平衡模拟，，，，检查是否有density.rst和density.mdcrd文件生成 ,命令如下
：

pmemd.cuda -O -i density.in -o density.out -p com_solvated.prmtop -c heat.rst -r density.rst -x density.mdcrd -ref heat.rst

所用参数文件 density.in,

heat ras-raf&cntrlimin=0,irest=1,ntx=5,nstlim=25000,dt=0.002,ntc=2,ntf=2,cut=8.0, ntb=2, ntp=1, taup=1.0,ntpr=500, ntwx=500,ntt=3, gamma_ln=2.0,temp0=300.0, ig=-1,ntr=1, restraintmask=':1-242',restraint_wt=2.0,/

• 体系平衡

最后进行体系平衡模拟，
，检查是否有equil.rst和equil.mdcrd生成

pmemd.cuda -O -i equil.in -o equil.out -p com_solvated.prmtop -c density.rst -r equil.rst -x equil.mdcrd

equil.in参数文件内容如下：

heat ras-raf&cntrlimin=0,irest=1,ntx=5,nstlim=250000,dt=0.002,ntc=2,ntf=2,cut=8.0, ntb=2, ntp=1, taup=2.0,ntpr=1000, ntwx=1000,ntt=3, gamma_ln=2.0,temp0=300.0, ig=-1,/

利用process_mdout.pl脚本检查上述平衡是否达到合理水平，，命令如下

process_mdout.pl heat.out density.out equil.out

利用xmgrace作图：

xmgrace summary.DENSITYxmgrace summary.TEMPxmgrace summary.ETOT

温度平衡结果

密度平衡结果

体系平衡结果

• 成品模拟
  

从上面分析可知体系基本上达到平衡，，，因此可以进行成品模拟步骤，，，，命令如下

pmemd.cuda -O -i prod.in -o prod.out -p com_solvated.prmtop -c equil.rst -r prod1.rst -x prod1.mdcrd

prod.in参数为
：

heat ras-raf&cntrlimin=0,irest=1,ntx=5,nstlim=500000,dt=0.002,ntc=2,ntf=2,cut=8.0, ntb=2, ntp=1, taup=2.0,ntpr=5000, ntwx=5000,ntt=3, gamma_ln=2.0,temp0=300.0, ig=-1, /

同样的利用Chimera的MD movie工具打开模拟结果文件的prod1.mdcrd和com_solvated.prmtop，，观看动力学模拟动画效果。
。

Amber动力学模拟动态展示效果

通过MD movie--Analysis--Plot--RMSD工具绘制体系的RMSD图,分析结果表明整个体系在模拟的时间段内处于2.0 Å范围内波动，
，，，是合理范围，，可用于后续结合自由能计算。。。。

动力学模拟体系RMSD结果

Amber集成了很多种结合自由能计算，，，例如MMGBSA或者MMPBSA,这里我们用MMGBSA来计算Abl和伊马替尼的结合自由能，，Amber中计算结合自由能的模块为MMPBSA.py命令如下：

python $AMBERHOME/bin/MMPBSA.py -O -i mmpbsa.in -o result.dat -sp com_solvated.prmtop -cp com.prmtop -rp receptor.prmtop -lp ligand.prmtop -y prod*.mdcrd

mmpbsa.in参数文件指定了计算结合自由能类型，，具体参数内容如下：

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

运行完毕会生成一系列记载能量信息的文件，，，最终结果包含在命令指定的dat后缀的文件中，，可提取数值进行作图分析
。。。。通过计算可知Abl蛋白质与imatinib的结合自由能为-62.73 kcal/mol,具体能量项如下图所示，
，主要包括了范德华相互作用、、静电相互作用、、溶剂自由能等
。。。。

Abl和imatinib结合自由能计算结果

上述教程是在舜王生物云计算平台上运行的细节，
，，在计算平台上提交作业
，，请不要在管理节点逐条需要消耗计算资源的命令，，管理节点只有2核4G内存的计算资源，，进行动力学模拟必须将作业提交到运算节点进行。。。。

因此需要将上述命令写进脚本中，，然后利用如下命令提交即可：

sbatch -N 1 -p g-v100-1 -c 12 amber.sh

N为节点的数量，
，
，，这里输入的是1。
。。
。p为选择的PARTITION，，这里使用的是V100卡（g-v100-1）。。。。查看作业运行情况及参数详细介绍可使用slurm命令进行查看。。

amber.sh内容是上述教程中所有命令的脚本形式，，
，，内容为：

#!/bin/bashmodule add Anaconda3/2020.02source /public/software/.local/easybuild/software/amber/aber20/amber.shulimit -s unlimitedulimit -l unlimitedantechamber -i ligand.mol2 -fi mol2 -o ligand.ante.mol2 -fo mol2parmchk2 -i ligand.ante.mol2 -f mol2 -o ligand.frcmodpdb4amber -i protein.pdb -o rec.pdb --reduce --dry -y --nohydtleap -s -f tleap.inpmemd.cuda -O -i min.in -o min.out -p com_solvated.prmtop -c com_solvated.inpcrd -r min.rst -ref com_solvated.inpcrdpmemd.cuda -O -i heat.in -o heat.out -p com_solvated.prmtop -c min.rst -r heat.rst -x heat.mdcrd -ref min.rstpmemd.cuda -O -i density.in -o density.out -p com_solvated.prmtop -c heat.rst -r density.rst -x density.mdcrd -ref heat.rstpmemd.cuda -O -i equil.in -o equil.out -p com_solvated.prmtop -c density.rst -r equil.rst -x equil.mdcrdpmemd.cuda -O -i prod.in -o prod1.out -p com_solvated.prmtop -c equil.rst -r prod1.rst -x prod1.mdcrdpython $AMBERHOME/bin/MMPBSA.py -O -i mmpbsa.in -o result.dat -sp com_solvated.prmtop -cp com.prmtop -rp receptor.prmtop -lp ligand.prmtop -y prod*.mdcrd

上述脚本中前面5行命令是告诉舜王生物云计算平台需要调用Amber程序，，，，后面的命令就是本教程所使用的具体命令了。。。

因此只需要将protein.pdb和ligand.mol2文件、、、所有in参数文件上传到舜王生物云计算平台，，，，然后提交amber.sh作业脚本即可完成本教程所有内容。
。。本教程参数文件获取地址：

 https://pan.baidu.com/s/1n3HuK9dp9o_uTdBjLuPpdQ   提取码: 49x8

• 高性能计算资源，，极大的节省计算成本

• 支持海量的计算工具，，，，且开箱即用

• 计算任务进度实时追踪提示的贴心服务

•  详细耐心的技术咨询