[MCPB.py踩坑]当金属蛋白中参与配位的残基数量>9时需要手动修改残基名

fffff

我在参考基于amber的MCPB.py构建含金属离子蛋白的力场参数 - 分子模拟 (Molecular Modeling) - 计算化学公社 (keinsci.com)使用MCPB.py构建金属蛋白拓扑文件用于分子动力学模拟。目前量化部分已经全部计算完毕，已生成tleap输入文件(见附件)。在执tleap -s -f 1DFS_tleap.in > 1DFS_tleap.out生成拓扑时遇到如下报错

1. /home/xxx/envs/AmberTools23/bin/teLeap: Fatal Error!
   
2. bond: Argument #2 is of type String must be of type: [atom]
   
3.     Here are some suggestions for correcting this error:
   
4.     Verify the type of an argument with the desc command.
   
5.     Check for alternate argument names with the list command.
   
6. usage:  bond <atom1> <atom2> [order]
   

复制代码

似乎是bond项的残基号和原子名有问题，首先检查载入的mol2信息

1. CM1 = loadmol2 CM1.mol2
   
2. CM2 = loadmol2 CM2.mol2
   
3. CM3 = loadmol2 CM3.mol2
   
4. CM4 = loadmol2 CM4.mol2
   
5. CM5 = loadmol2 CM5.mol2
   
6. CM6 = loadmol2 CM6.mol2
   
7. CM7 = loadmol2 CM7.mol2
   
8. CM8 = loadmol2 CM8.mol2
   
9. CM9 = loadmol2 CM9.mol2
   
10. CM10 = loadmol2 CM10.mol2
    
11. CM11 = loadmol2 CM11.mol2
    
12. CD1 = loadmol2 CD1.mol2
    
13. CD2 = loadmol2 CD2.mol2
    
14. CD3 = loadmol2 CD3.mol2
    
15. CD4 = loadmol2 CD4.mol2

复制代码

11个CYM和4个CD，无误。
然后，检查bond命令内容和通过loadpdb载入的pdb文件(1DFS_mcpbpy.pdb)中的残基号和原子名对应关系，参与配位的S负离子原子名为SG，镉离子的原子名为CD。四个镉离子的残基号是32-35，也都能对得上(pdb文件和in文件的残基号和原子名对得上）

1. bond mol.3.SG mol.33.CD
   
2. bond mol.4.SG mol.33.CD
   
3. bond mol.4.SG mol.35.CD
   
4. bond mol.6.SG mol.35.CD
   
5. bond mol.7.SG mol.34.CD
   
6. bond mol.7.SG mol.35.CD
   
7. bond mol.11.SG mol.34.CD
   
8. bond mol.14.SG mol.33.CD
   
9. bond mol.14.SG mol.34.CD
   
10. bond mol.18.SG mol.33.CD
    
11. bond mol.20.SG mol.32.CD
    
12. bond mol.20.SG mol.35.CD
    
13. bond mol.27.SG mol.32.CD
    
14. bond mol.29.SG mol.32.CD
    
15. bond mol.30.SG mol.32.CD
    
16. bond mol.30.SG mol.34.CD

复制代码

然后通过tleap的desc命令检查，发现了端倪，即tleap读入的和我直接打开pdb看到的不一致

从29号开始，残基名就没有被正确识别。此处CM10代表CYM10（体系中的第10个CYM），CM11代表CYM11（体系中的第11个CYM），本来29号是CM10,30号是CM11，但是此处把29号识别为了CM1（CYM1体系中第一个CYM）30号CM11没有被识别，C末端的ALA被放到了30位，导致CD的序号都少1。

原因在于pdb文件中残基名一般是三字符，后加空格再加链标识符，即第18-20列是残基名，第21列空格，第22列链标识符。但是MCPB.py在生成xxx_mcpbpy.pdb文件时，对于四字符的CM10和CM11，仍然按照[残基名][空格][链标识符]处理，导致链标识符从22列移动到了23列，导致后面的残基号也不在正确的列，进而导致无法被amber的tleap正确识别。

于是手动删除多余空列，让链标识符回到第22列，在telap中执行desc命令仍有问题

此时残基号可以被正确识别，29和30号是CYM，31是ALA，32开始是CD。但此时如果直接继续执行tleap生成拓扑仍会报错，提示缺参数。出现该问题是因为tleap在读取pdb的残基名时，只读前三个字符，即CM10只读取CM1，CM11只读取CM1，和真正的第一个CYM CM1重复，导致参数缺失。解决措施是把CM10改成CMX，把CM11改成CMY，使用tleap的desc命令可以看到残基名和残基号均被正确识别

此外，还要修改前面生成CYM残基mol2文件相应位置，将其CM10修改为CMX，CM11修改为CMY，同时将telap.in文件中load字段做相应修改，tleap才可以正常识别。

1. @<TRIPOS>MOLECULE
   
2. CMX
   
3.    10     9     1     0     0
   
4. SMALL
   
5. RESP Charge
   
6. 
   
7. 
   
8. @<TRIPOS>ATOM
   
9.       1 N          -7.3020   -4.0340    2.9680 N         1 CM10    -0.415700
   
10.       2 H          -6.8280   -3.1430    2.9760 H         1 CM10     0.120284
    
11.       3 CA         -6.3850   -5.3090    2.5050 CX        1 CM10    -0.035100
    
12.       4 HA         -7.0000   -6.1990    2.6420 H1        1 CM10     0.075008
    
13.       5 CB         -6.0920   -5.1860    0.8710 CT        1 CM10    -0.015327
    
14.       6 HB2        -5.5630   -6.0920    0.5730 H1        1 CM10     0.088956
    
15.       7 HB3        -7.0570   -5.2020    0.3650 H1        1 CM10     0.088956
    
16.       8 SG         -5.1260   -3.7270    0.1890 Y3        1 CM10    -0.587176
    
17.       9 C          -5.0040   -5.5720    3.3940 C         1 CM10     0.597300
    
18.      10 O          -4.5500   -6.7520    3.4100 O         1 CM10    -0.567900
    
19. @<TRIPOS>BOND
    
20.      1    1    2 1
    
21.      2    1    3 1
    
22.      3    3    4 1
    
23.      4    3    5 1
    
24.      5    3    9 1
    
25.      6    5    6 1
    
26.      7    5    7 1
    
27.      8    5    8 1
    
28.      9    9   10 1
    
29. @<TRIPOS>SUBSTRUCTURE
    
30.      1 CMX         1 TEMP              0 ****  ****    0 ROOT
    

复制代码

tleap.in文件

1. CM1 = loadmol2 CM1.mol2
   
2. CM2 = loadmol2 CM2.mol2
   
3. CM3 = loadmol2 CM3.mol2
   
4. CM4 = loadmol2 CM4.mol2
   
5. CM5 = loadmol2 CM5.mol2
   
6. CM6 = loadmol2 CM6.mol2
   
7. CM7 = loadmol2 CM7.mol2
   
8. CM8 = loadmol2 CM8.mol2
   
9. CM9 = loadmol2 CM9.mol2
   
10. CMX = loadmol2 CMX.mol2
    
11. CMY = loadmol2 CMY.mol2
    
12. CD1 = loadmol2 CD1.mol2
    
13. CD2 = loadmol2 CD2.mol2
    
14. CD3 = loadmol2 CD3.mol2
    
15. CD4 = loadmol2 CD4.mol2

复制代码

后续可继续按照教程检查拓扑，转换gmx拓扑等。