Back

gromacs中构建非标准氨基酸残基力场

Gromacs非标准残基构建

李老师其实已经有非常详细的实例来讲解如何进行非标准残基的搭建方法,此方法是之前根据李老师的方法结合amber教程进行改编的方法,可以作为参考。

我们使用ser磷酸化的蛋白作为实验,由于还未发表文章,暂时不提供示例步骤下载。若有问题可以邮件联系我交流:kangsgo#vip.qq.com

1.单独保存非标准残基

我们首先使用pymol将sep独立导出

1
2
select not resn SEP
remove sele

将其保存为sep.pdb文件,当然也可以直接用gedit等打开文件,将SEP部分(HETAM)copy导出

sep文件内容如下:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
ATOM      1  N   SEP A 219      10.386  -7.407  38.687  1.00 31.68           N  
ATOM      2  CA  SEP A 219      10.979  -8.567  39.325  1.00 31.58           C  
ATOM      3  C   SEP A 219      12.450  -8.306  39.481  1.00 28.97           C  
ATOM      4  O   SEP A 219      13.189  -8.279  38.504  1.00 28.84           O  
ATOM      5  CB  SEP A 219      10.724  -9.808  38.495  1.00  0.00           C  
ATOM      6  OG  SEP A 219      11.474 -10.889  39.010  1.00  0.00           O  
ATOM      7  P   SEP A 219      11.419 -12.337  38.313  1.00  0.00           P1+
ATOM      8  O2P SEP A 219      12.339 -13.139  39.147  1.00  0.00           O1-
ATOM      9  O1P SEP A 219       9.990 -12.703  38.422  1.00  0.00           O1-
ATOM     10  O3P SEP A 219      11.893 -12.061  36.938  1.00  0.00           O1-
ATOM     11  H   SEP A 219      10.767  -7.078  37.811  1.00 31.68           H  
ATOM     12  HA  SEP A 219      10.534  -8.738  40.305  1.00 31.58           H  
ATOM     13  HB2 SEP A 219       9.662 -10.058  38.502  1.00  0.00           H  
ATOM     14  HB3 SEP A 219      11.035  -9.614  37.467  1.00  0.00           H  

不难看出带的电荷为-2

我们需要补全两端的氢键(封端),我是用gview补全的,也可以用pymol,chimera等补全,但是需要注意的是要查看补全的氢键是否正确,我在chimera上补全的氢键就有问题,对于补全氢键,个人的经验觉得正确度排序为:gv>chimera>pymol

补全后如下:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
HETATM    1  N   SEP A 219      10.386  -7.407  38.687                       N
HETATM    2  CA  SEP A 219      10.979  -8.567  39.325                       C
HETATM    3  C   SEP A 219      12.450  -8.306  39.481                       C
HETATM    4  O   SEP A 219      13.189  -8.279  38.504                       O
HETATM    5  CB  SEP A 219      10.724  -9.808  38.495                       C
HETATM    6  OG  SEP A 219      11.474 -10.889  39.010                       O
HETATM    7  P   SEP A 219      11.419 -12.337  38.313                       P
HETATM    8  O2P SEP A 219      12.339 -13.139  39.147                       O
HETATM    9  O1P SEP A 219       9.990 -12.703  38.422                       O
HETATM   10  O3P SEP A 219      11.893 -12.061  36.938                       O
HETATM   11  H   SEP A 219      10.767  -7.078  37.811                       H
HETATM   12  HA  SEP A 219      10.534  -8.738  40.305                       H
HETATM   13  HB2 SEP A 219       9.662 -10.058  38.502                       H
HETATM   14  HB3 SEP A 219      11.035  -9.614  37.467                       H
HETATM   15 2H   SEP A 219      10.544  -6.657  39.329                       H
HETATM   16  H   SEP A 219      12.857  -8.142  40.457                       H
TER      16      SEP A 219
END
CONECT    1    2   11   15
CONECT    2    1    3    5   12
CONECT    3    2    4   16
CONECT    4    3
CONECT    5    2    6   13   14
CONECT    6    5    7
CONECT    7    6    8    9   10
CONECT    8    7
CONECT    9    7
CONECT   10    7
CONECT   11    1
CONECT   12    2
CONECT   13    5
CONECT   14    5
CONECT   15    1
CONECT   16    3

2.制作高斯输入文件

1
antechamber -fi pdb -i sep.pdb -fo gcrt -o sep.gjf -ch "sep.chk" -gm "%mem=2048MB" -gn "%nproc=4" -nc -2

如果按照李老师的方法还需要使用-ge 生成高斯esp文件,通过iop(6/50=1),可以参考下面一条命令

1
antechamber -fi pdb -i sep.pdb  -o ligand.gjf -fo gcrt -pf y -gn "%nproc=8" -gm "%mem=1000MB" -ch "ligand" -gk "#HF/6-31G* SCF=tight Test Pop=MK iop(6/33=2) iop(6/42=6) iop(6/50=1) opt" -ge ligand.gesp -gv 1

命令解释如下: -fi为输入格式 -i为输入文件 -fo为输出文件格式 gcrt为高斯gjf格式 -ch为高斯chk文件 -gm -gn为使用内存和核数 -nc 为电荷数

3.转化为gromacs文件

实际上与蛋白配体中的小分子制作方式是一样的,首先在这里下载ACPYPE,注意需要下载sf版本,否则而面角的表达方式不对,首先生成mol2文件

1
antechamber -fi gout -i sep.log -rn SEP -fo mol2 -o sep.mol2 -c resp -pf y -at amber

在这里最好修改一下原子名称,改成与sep.pdb对应,且需要对其进行修改,因为gromacs加氢从1开始,故将HB2,HB3改为HB1,HB2 最后的结果如下:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
...
@<TRIPOS>ATOM
      1 N            3.0030    -1.1500     0.1620 N         1 SEP      -1.141452
      2 CA           1.8370    -0.4940    -0.4350 CT         1 SEP       0.462377
      3 C            1.8020     0.9980    -0.1460 C          1 SEP       0.457721
      4 O            2.8040     1.5940     0.1700 O          1 SEP      -0.633121
      5 CB           0.5190    -1.1870    -0.0550 CT         1 SEP       0.192426
      6 OG          -0.5100    -0.6800    -0.7990 OS         1 SEP      -0.641051
      7 P           -1.7850     0.1390     0.0880 P          1 SEP       1.339193
      8 O2P         -2.7260     0.5020    -1.0190 O2         1 SEP      -0.920441
      9 O1P         -2.2370    -0.9110     1.0650 O2         1 SEP      -0.920441
     10 O3P         -1.0110     1.2850     0.7110 O2         1 SEP      -0.920441
     11 H            2.7710    -1.3680     1.1160 H          1 SEP       0.367588
     12 HA           1.9150    -0.5690    -1.5210 H1         1 SEP      -0.041482
     13 HB1          0.6570    -2.2570    -0.2470 H1         1 SEP      -0.007243
     14 HB2          0.3490    -1.0790     1.0140 H1         1 SEP      -0.007243
     15 H2           3.7290    -0.4600     0.2230 H          1 SEP       0.367588
     16 H1           0.8340     1.4820    -0.2400 HA         1 SEP       0.046022
....

生成的mol2文件使用parmchk2进行补全参数

1
parmchk2 -i sep.mol2 -f mol2 -o sep.frcmod

创建leap.in文件,输入如下内容(注意sep修改成自己的名称):

1
2
3
4
5
6
7
source leaprc.ff99SBildn
source leaprc.gaff
loadamberparams sep.frcmod
lig=loadmol2 sep.mol2
check lig
saveamberparm lig sep.prmtop sep.inpcrd
quit

运行

1
tleap -f leap.in

运行acpype

1
python2.7 acpype.py -p sep.prmtop -x sep.inpcrd -d

将会生成SEP_GMX.topSEP_GMX.gro两个文件,其中SEP_GMX.gro我们用不着

4.整理残基的rtp条目

为了保存力场rtp信息,我们需要创建一个sep.rtp文件 我们需要对SEP_GMX.gro文件进行整理,即把虚拟的原子电荷附加到N端头和C端尾部的位置。 如:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
[ atoms ]
;   nr  type  resi  res  atom  cgnr     charge      mass       ; qtot   bond_type
     1    N     1   SEP     N    1    -1.141453     14.01000 ; qtot -1.141
     2   CT     1   SEP    CA    2     0.462377     12.01000 ; qtot -0.679
     3    C     1   SEP     C    3     0.457721     12.01000 ; qtot -0.221
     4    O     1   SEP     O    4    -0.633122     16.00000 ; qtot -0.854
     5   CT     1   SEP    CB    5     0.192426     12.01000 ; qtot -0.662
     6   OS     1   SEP    OG    6    -0.641052     16.00000 ; qtot -1.303
     7    P     1   SEP     P    7     1.339196     30.97000 ; qtot 0.036
     8   O2     1   SEP   O2P    8    -0.920442     16.00000 ; qtot -0.884
     9   O2     1   SEP   O1P    9    -0.920442     16.00000 ; qtot -1.805
    10   O2     1   SEP   O3P   10    -0.920442     16.00000 ; qtot -2.725
    11    H     1   SEP     H   11     0.367588      1.00800 ; qtot -2.358
    12   H1     1   SEP    HA   12    -0.041482      1.00800 ; qtot -2.399
    13   H1     1   SEP   HB1   13    -0.007243      1.00800 ; qtot -2.406
    14   H1     1   SEP   HB2   14    -0.007243      1.00800 ; qtot -2.414
    15    H     1   SEP    H2   15     0.367588      1.00800 ; qtot -2.046
    16   HA     1   SEP    H1   16     0.046022      1.00800 ; qtot -2.000

1 N的电荷应该变为-1.141453+0.367588 3 C的电荷应该变为0.457721+0.046022

然后用我写的脚本进行内容的提取: 首先创建config.txt文件,在其中写入需要删除的虚拟原子(每一行一个)

1
2
H
HA

然后运行

1
python rtp.py SEP_GMX.top SEP>data.txt

结果如下:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
[ bondedtypes ]
; bonds  angles  dihedrals  impropers all_dihedrals nrexcl HH14 RemoveDih
     1       1          9          4        1         3      1     0
[ SEP ]
 [ atoms ]
      N     N         -0.773865     1
...
 [ bonds ]
      N    CA        1.4620e-01     2.7506e+05
      N     H        1.0130e-01     3.3740e+05
     CA     C        1.5240e-01     2.6192e+05
...
 [ angles ]
      N    CA     C      1.0906e+02        5.6066e+02
      N    CA    CB      1.1161e+02        5.5153e+02
      N    CA    HA      1.0888e+02        4.1706e+02
     CA     N     H      1.1768e+02        3.8325e+02
     CA     C     O      1.2320e+02        5.6400e+02
...
 [ dihedrals ] ; propers
     N    CA     C         O    180.00   0.00000     9
     N    CA    CB        OG      0.00   0.65084     9
     N    CA    CB       HB1      0.00   0.65084     9
     N    CA    CB       HB2      0.00   0.65084     9
...
 [ impropers ]

我们需要在[bonds]中增加前一个氨基酸残基的联系连接,[impropers]中增加前后氨基酸残基的联系连接,可以简单的看非标准残基来源的氨基酸残基模板(aminoacids.rtp)或者查看SEP_GMX.top文件中删除的两个氢原子的部分,最后整理结果如下:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
...
 [ bonds ]
     -C     N
      N    CA        1.4620e-01     2.7506e+05
      N     H        1.0130e-01     3.3740e+05
     CA     C        1.5240e-01     2.6192e+05
     CA    CB        1.5380e-01     2.5179e+05
...
 [ impropers ]
    -C    CA     N     H
    CA    +N     C     O

5.整理残基的hdb条目

由于氢键一般有问题,我们需要重新绘制氢键,那么我们创建一个sep.hdb的文件,按照下图来进行设置(图片来源李老师博客): 图

结果如下图(可以抄SER的信息):

1
2
3
4
SEP  3
1   1   H    N  -C   CA
1   5   HA   CA  N   CB    C
2	6	HB	CB	CA	OG	

6.模拟尝试

将gromacs 目录下(gromacs/share/gromacs/top)中的amber99sb-ildn.ff拷贝至实验目录下(如果找不到可以用which gmx_mpi),将sep.rtp文件和sep.hdb文件放到拷贝出来的amber99sb-ildn.ff目录下。且在top目录下的residuetypes.dat中申明SEP为蛋白,即:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
SEP     Protein ; sep,增加部分
ABU	Protein
ACE	Protein
AIB	Protein
ALA  	Protein
ARG  	Protein
ARGN	Protein
ASN	Protein
ASN1  	Protein
...

最后进行模拟尝试看是否可以跑通

1
gmx_mpi pdb2gmx -f model.pdb -o model_process.gro -water tip3p -ignh

若提示前一个残基或者末端残基报错,那么检查是否是上述两个文件出错,若提示缺失(SEP的)原子,那么创建sep.atp文件,在SEP_GMX.top,找到提示的缺失原子,放入该文件中,最后将sep.atp文件放到拷贝出来的amber99sb-ildn.ff目录下。

参考资料: GROMACS非标准残基教程2-芋螺毒素小肽实例 向gromacs中添加小分子力场方法 amber中非标准氨基酸残基的参数生成 Simulating the Green Fluorescent Protein 使用AmberTools+ACPYPE+Gaussian创建小分子GAFF力场的拓扑文件

Licensed under CC BY-NC-SA 4.0