软件与材料
本文内容使用到如下软件:
- Ambertools 20
- GROMACS 2019.5
- Gaussian 16
如果使用win10操作系统,可以自己编译前两个软件,也可以找别人编译好的版本。比如Ambertools、GROMACS2019.5;在此表示感谢。Gaussian 16 是要付费购买的,在本文中是用于计算配体的RESP电荷,RESP电荷还可以通过Multiwfn得到;配体还可以采用bcc电荷(准确性不如RESP),可以直接用Ambertools计算得到。
本文还用到一个蛋白protein.pdb和一个从drugbank下载的配体小分子ligand.pdb。蛋白将利用GROMACS的力场处理,而配体将用amber的GAFF力场处理。
配体处理
配体加氢
本次下载得到的ligand.pdb文件里是没有包含非极性氢的,所以需要给它加氢。加氢的方式有很多,比如Ambertools里的reduce、或者pymol等。
reduce ligand.pdb > ligand_h.pdb
reduce 可能更适用于氨基酸组成的分子;不论使用什么软件进行加氢操作,最后一定要检查结果文件中是否包含了加上的氢。还需要检查加氢之后的结构是否合理,不合理的结构无法通过高斯进行计算。
之后通过Ambertools的antechamber处理成Gaussian的输入文件。
antechamber -i ligand_h.pdb -fi pdb -o ligand.gjf -fo gcrt -pf y
# 当然输入的配体文件也可以是mol2等
# antechamber -i liand_h.mol2 -fi mol2 -o ligand.gjf -fo gcrt -pf y
然后利用Gaussian 16 计算:
g16 < ligand.gjf > ligand.out
在Gaussian的输出文件中就有ESP电荷数据;利用antechamber得到包含RESP电荷的mol2文件:
antechamber -i ligand.out -fi gout -c resp -o ligand_resp.mol2 -fo mol2 -pf y
要在Ambertools中对配体应用GAFF力场,我们还需要生成一个输入文件:
parmchk -i ligand_resp.mol2 -f mol2 -o ligand.frcmod
ligand.frcmod中包含小分子的一些必要参数。
如此,我们所需要的的Ambertools输入文件就准备好了。
对配体应用GAFF力场
打开Ambertools中的tleap,输入如下命令:
# 载入力场
source leaprc.gaff
# 载入配体参数化文件
loadamberparams ligand.frcmod
# 载入配体
lig = loadmol2 ligand_resp.mol2
# 检查配体
check lig
# 保存amber输入文件
saveamberparm lig lig.prmtop lig.inpcrd
转换配体结构和拓扑文件到GROMACS格式
接下来我们需要将amber的输入文件转换为GROMACS可用的输入文件,也即top文件和gro文件。
通过Ambertools中的acpype程序进行转换:
# 输出帮助信息
# acpype -h
# 做文件转换,输出文件名默认为Ambertools中unit的名字
acpype -p lig.prmtop -x lig.inpcrd
如此我们就得到了top文件和gro文件;配体的top文件还需要一点修改才可以被GROMACS使用。
构建模拟体系
生成蛋白质拓扑文件
我们已经有了蛋白文件protein.pdb,配体的top和gro文件;接下来首先要做的就是创建蛋白质的拓扑文件:
gmx pdb2gmx -f protein.pdb -o protein.gro -water spc -ignh
# 选择力场
# 这里选择amber99sb力场
然后我们在相应目录下就可以看到生成了protein.gro、topol.top、posre.itp。topol.top里面包含了选择的力场的一些参数以及蛋白质的拓扑结构参数等等。
接下来我们需要做一系列文件更改,来把配体加到我们的体系中。
拓扑文件修改
先说明本例子中各文件的修改步骤:
- 复制protein.gro并重命名为complex.gro
- 将配体的gro文件的坐标部分,复制到complex.gro的坐标部分后面,并相应修改第二行的原子数目
- 将配体的top文件重命名为ligand.itp
- 将pdb2gmx命令生成的topol.top文件重命名为protein.itp
- 新建文件topol.top
- 将protein.itp中
[ moleculetype ]
部分之前的内容完全剪切到topol.top - 将protein.itp中
[ dihedrals ]
部分之后的内容完全剪切到topol.top,包括各种文件的导入部分、[ system ]以及[ molecules ]。 - 在topol.top中导入蛋白质位置限制文件的前面,导入蛋白质拓扑文件protein.itp
- 对照ligand.itp中最后的
[ molecules ]
,在topol.top中[ molecules ]
部分添加配体的名字和数目;在topol.top导入水拓扑文件的前面导入配体拓扑文件ligand.itp。 - 删去ligand.itp中最后的
[ system ]
和[ molecules ]
部分 - 适当修改topol.top中[ system ]的名字(这个可能无所谓)
- 对于ligand.itp中
[ moleculetype ]
前面的内容,我们也需要整合到topol.top中,一般包含两部分,[ defaults ]
和[ atomtypes ]
。对于这两部分内容,直接复制到topol.top文件中导入蛋白质拓扑文件部分的前面即可。 - 现在我们topol.top文件的最开始,是导入力场参数文件,之后就是配体的
[ defaults ]
部分,由于力场文件中通常也包含这一部分,所以会出错,将[ defaults ]
部分注释掉;然后尝试对整个体系进行溶剂化、电荷平衡操作,如果不报错,可能就可以了,如果报错,那可能还需要进一步修改topol.top的力场部分。比如用配体的[ defaults ]
部分置换掉力场的[ defaults ]
部分等。 - 有关力场文件的详细信息,可参考jerkwin的博文,主要是5.7.1部分;在此表示感谢。
关于本体系中拓扑文件修改的详细情况,说明如下。
首先,在配体和蛋白质的拓扑文件(.top)中,主要包含如下几个部分的内容:
- [ defaults ]
- [ atomtypes ]
- [ moleculetype ]
- [ atoms ]
- [ bonds ]
- [ pairs ]
- [ angles ]
- [ dihedrals ]
- [ system ]
- [ molecules ]
前两个部分[ defaults ]
和[ atomtypes ]
主要是相应的力场参数的部分,后面的部分则定义了一些分子拓扑结构相关的参数以及系统参数。一般的itp文件的内容则主要包括3-8的部分,所以我们把ligand.top重命名为ligand.itp、pdb2gmx生成的topol.top重命名为protein.itp之后需要去掉多余的部分;但是去掉的部分并不是不需要的,这部分内容需要与体系的topol.top文件进行整合。
首先将protein.itp的1-2、9-10部分剪切到空白的topol.top中,并include蛋白质的itp文件,得到topol.top如下:
; Include forcefield parameters
#include "amber99sb.ff/forcefield.itp"
; Include protein topology
#include "protein.itp"
; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif
; Include water topology
#include "amber99sb.ff/spc.itp"
#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
#endif
; Include topology for ions
#include "amber99sb.ff/ions.itp"
[ system ]
; Name
Protein
[ molecules ]
; Compound #mols
Protein 1
然后我们需要将ligand.itp中的内容整合到topol.top中来。
首先是整合[ system ]
和[ molecules ]
部分,整合后的topol.top的这两部分如下:
[ system ]
; Name
Protein_and_ligand
[ molecules ]
; Compound #mols
Protein 1
MOL 1
;主要就是添加了配体这一行,这里配体名为MOL
整合之后记得删除ligand.itp中的这两部分,并在topol.top导入水拓扑文件的前面导入ligand.itp。
然后整合ligand.itp的前两个部分到topol.top中,直接剪切并粘贴到topol.top导入蛋白质拓扑文件前面、导入力场的后面,然后注释掉[ default ]
部分。
整合之后的topol.top如下:
; Include forcefield parameters
#include "amber99sb.ff/forcefield.itp"
; [ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
; 1 2 yes 0.5 0.8333
[ atomtypes ]
;name bond_type mass charge ptype sigma epsilon Amb
c2 c2 0.00000 0.00000 A 3.39967e-01 3.59824e-01 ; 1.91 0.0860
c1 c1 0.00000 0.00000 A 3.39967e-01 8.78640e-01 ; 1.91 0.2100
os os 0.00000 0.00000 A 3.00001e-01 7.11280e-01 ; 1.68 0.1700
o o 0.00000 0.00000 A 2.95992e-01 8.78640e-01 ; 1.66 0.2100
; Include protein topology
#include "protein.itp"
; Include Position restraint file
#ifdef POSRES
#include "posre.itp"
#endif
; Include ligand topology
#include "ligand.itp"
; Include water topology
#include "amber99sb.ff/spc.itp"
#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
#endif
; Include topology for ions
#include "amber99sb.ff/ions.itp"
[ system ]
; Name
Protein_and_ligand
[ molecules ]
; Compound #mols
Protein 1
MOL 1
后续建模
通过GROMACS进行后续的建模,如果在添加离子等后续步骤中不报错的话,体系一般就是ok的。
gmx editconf -f complex.gro -o newbox.gro -bt cubic -d 1.0
gmx solvate -cp newbox.gro -cs spc216.gro -p topol.top -o solv.gro
gmx grompp -f em.mdp -c solv.gro -r solv.gro -p topol.top -o ions.tpr -maxwarn 2
gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -np 5
关于力场
这种体系中,小分子的力场最好选择和蛋白质的力场比较相近的力场;也即同一个体系中,所有的物质,最好都采用同一个力场或者同一类力场。本文中,GAFF力场和amber99sb算是同一类力场吧。
如果采用不同类型的力场(如假设本文中蛋白质采用GROMOS系列的力场),可能会有一些理论上的问题,请咨询相关的资深人员。
关于多配体体系
如果蛋白配体体系中需要有多个位置不同的配体,那可以通过pymol等工具,载入ligand_resp.mol2文件,通过复制移动等功能实现蛋白周围的不同位置配体的生成,然后分别保存每个配体为mol2文件。这样就可以得到多个有resp电荷的配体文件,但要注意:pymol等工具导出的文件中原子类型可能和ligand_resp.mol2中定义的「原子类型」不同了,需要对照ligand_resp.mol2文件对原子类型一列进行修改,才能与ligand.frcmod配合进行后续建模。之后再通过Ambertools处理,acpype转换。最后整合和include到体系中应该就可以了。
关注公众号蓝色火烈鸟,了解更多技术笔记和工具推荐: