对有非标准残基的蛋白质做蛋白质-配体分子动力学

lisanoid

      在21年12月份参加了sob社长的分子动力学与gromacs程序培训班，受益匪浅。之后回来了自己做一些工作，主要是围绕小分子化合物和蛋白质结合（分子对接）和分子动力学模拟（MD）开展的。    在做了大概几十个蛋白-配体复合物的MD之后。总体上感觉，以带有非标准残基的蛋白质构建工作最繁琐，正好这次的蛋白质又碰上该问题，就借公社宝地一用，记录一次处理全过程。如果烂尾了，就请版主把此帖删掉。
    首先介绍一下问题背景：在使用gromacs做蛋白质的动力学模拟时，主要使用pdb2gmx命令来构建跑MD所需的拓扑文件，但当蛋白质中存在pdb2gmx不能识别的残基时（能识别的残基可以参见gromacs文件夹/share/gromacs/top/residuetypes.dat 中内容），则出现报错。解决该问题的中文教程，目前比较全面的有李继存老师博文，以及知乎上天津大学supernova的Fmoc残基生成教程。以及社长在论坛中发布的部分帖子。记录不太详细，工作繁琐，而且经常报错，没有什么简单可靠的方法。
   根据之前的非标准残基处理经验，主要内容参考http://sobereva.com/soft/Sobtop/   中的FAQ5 ，我把这块工作分了几步，分别是：（1） 裁剪残基，封装后形成新分子，做优化和振动分析；（2）用波函数分析软件multiwfn计算resp电荷（http://sobereva.com/441） ；（3）用multiwfn创建残基的rtp文件 （4）根据残基结构裁剪rtp文件 （5）把文件整合进gromacs相应文件中，在gromacs中做pdb2gmx处理生成蛋白质拓扑文件；（6）报错及修复过程 （7）MD后续及结果分析
   本帖拟分步记录处理过程，由于本人水平较低，处理过程出现了问题，正好请社长和多多指教。

BobXu

我也手搓了一个赖氨酸乳酰化的rtp和hdb文件，发现了一个大问题，就是在生成的rtp文件中，由于氢原子太多，我把HA， HB，HG，HD，HZ用到了赖氨酸上，但是乳酸上的氢原子不知道如何命名。后来发现生成的rtp文件中，氢原子全部是以H1，H2这样的数字标识的，我就按照赖氨酸的命名方式把rpt文件中的氢原子全部手动修改为HB1，HB2，HG1，HG2这种形式，在乳酸中的氢原子，我又重新定义了HL和HM来分别代表如乳酰羟基碳上的氢和甲基碳上的氢~~~~问了好多，都没回复，手动试一下，没想到成功了。还有一个，我是试了CHARM-GUI生成的文件，但是它生成的不对，无法跑MD，而且CHARM-GUI的乳酸化到最后都变成了LYN，并且报错。

Sam_

lisanoid 发表于 2022-6-27 20:25
分别是不含磷酸链nophos，含磷酸肽链withphos，以及配体分子的结构文件

删除磷酸链以后的蛋白肽链就没有这两个非标准残基了吧

dayu8278

大佬好，你提到的李继存老师博文能不能分享下

clock1993

hongyi166_ 发表于 2022-9-11 19:48
大佬讲的很全面，我基本上同样的问题也都遇到了一遍，只不过我的加入rtp之后出现HIS拓扑文件的缺失，尝 ...

大佬您好，请问您的这个问题解决了吗？最后是怎么处理的？谢谢！

znjnancy

lisanoid 发表于 2022-5-27 15:19
这次研究的目标蛋白是PRKC  蛋白质pdbID：3IW4。从RCSB网站（https://www.rcsb.org/）下载结构，
  采 ...

请问你这个gaussian的L1报错最后如何解决的呢

JohnCase

lisanoid 发表于 2022-5-29 18:40
3.构建rtp文件
    根据sobtop的FAQ5 ，rtp文件生成方法为：然后在Sobtop里载入模型分子的mol2文件，从chg ...

老师您好，请教一下这里选择了“4 Same as option 3, but missing ones are arbitrarily guessed”的话，就不需要让Sobtop基于Hessian算吧？

liu840397562

lisanoid 发表于 2022-5-29 19:39
4 rtp文件的编辑
rtp文件生成后，不能直接使用，因为当时生成时是加了甲胺基和乙酰基的，而正常的肽链中， ...

请问这步只能手动删除这些端基信息吗？ gromacs有无内置功能处理呢，我在网上搜索后没找到特别的答案，特来请教以下！

hongyi166_

大佬讲的很全面，我基本上同样的问题也都遇到了一遍，只不过我的加入rtp之后出现HIS拓扑文件的缺失，尝试了一下，发现没有磷酸化的同样蛋白不会出现这个报错，所以又重新调整了rtp文件，主要是改了原本的乙酰基和甲胺基为“-C和+N“，也就是和上下氨基酸残基的联系，然后报错就变成了”需要加入氢键相关部分，也就是ignh对应的需要gmx自己补氢“，再自己编写了hdb文件之后pdb2gmx这部就没有报错了。
现在的报错问题是，当时用乙酰和甲胺基做了两端的封堵后 算的resp电荷，现在删掉之后，残基的电荷不是整数了，请问楼主遇到了这个问题嘛？怎么解决呢？

lisanoid

分别是非标准氨基酸的rtp文件和hdb文件.

lisanoid

分别是不含磷酸链nophos，含磷酸肽链withphos，以及配体分子的结构文件

lisanoid

10. 任务的结束和总结
10.1  针对9中的报错处理，与社长请教后，发现是力场问题：使用sobtop处理得到的拓扑文件，不支持gromos力场，支持amber力场，而我生成残基之后选择了gromacs96 54a7力场，所以处理非标准氨基酸时所有的力场计算都出现了错误。把两个非标准氨基酸的数据拷到amber99sbildn力场里，问题解决。
10.2 继续计算时，出现新报错“Incorrect number of parameters - found 4, expected 3 or 6 for Proper Dih.  (after the function type).”
发现是sobtop生成拓扑文件后，再用pdb2gmx生成蛋白质topol文件，会出现二面角functype出现两次的问题。在topol文件中，找到对应提示行，把多出来的列删除即可，这个后面会出现很多次，所以可以把top文件，从开始提示的那一行开始，一直排查到最后。
10.3 能量最小化时，出现特定原子能量和受力异常大的情况，根据提示找到相应原子，在分子可视化软件里观察，发现pdb文件在经pulchra处理后，420号氨基酸val有氢原子与邻近的H原子重叠，这显然是不可接受的，重新优化蛋白质结构。
10.4 又出现受力异常情况，实在不胜其烦了，观察蛋白质结构以后，把肽链外围一条单独的链给剪去了，因为外围这条链问题最多，同时距离配体位置也比较远，干脆一删了之。
但理论上说，这是个磷酸化酶，剪去的短肽链中，有几个磷酸化位点，虽然距离比较远，但是磷酸化基团电性强，是不是有可能影响到长时间的动力学模拟呢？
目前水平有限，只能做到此种程度，现在把几个文件贴出来，请行家们指教，多谢！

lisanoid

9 报错处理（2）
昨天的No default Proper Dih. types 我想了一下，是不是之前构建top文件时，没有对非标准氨基酸加氢造成的？
于是对两个残基做了加氢操作，生成了hbd文件，并根据http://bbs.keinsci.com/thread-24334-1-1.html对H原子命名做了改进，。
重新生成蛋白质top文件后，仍然出现上次同样的错误，4500多个error....看来单纯加氢不行，可能还是命名或者编号出了某些问题导致的报错？
还有就是sobtop中生成的rtp文件中，有这么一句“; NOTICE: After generation of .top file using pdb2gmx, functype in [ dihedrals ]
;         will present twice, please manually keep the one given in this .rtp file”，这个要好好请教请教才行了。

lisanoid

8. 报错及处理（1）
上述文件汇总后，利用gromacs进行分子动力学处理，处理时遇到bug，具体操作是，当给复合物加水后，运行gmx grompp -f em.mdp -c complex_SOL.gro -p topol.top -o em.tpr -maxwarn 3命令时
报错：No default Proper Dih. types   共出现了4900多次
初步判断是非标准残基处理之后，部分信息被删掉，可能造成二面角信息不全。而且知乎上，天津大学supernova的Fmoc残基也出现了类似错误。因此把原来生成的TPO.rtp和SEP.rtp文件，重新补充信息后处理，但效果不彰，报错问题依然没有解决...而且报错中，很多都是标准氨基酸的4个原子之间的，也报错说没有二面角信息，因此怀疑是不是之前补充SEP和TPO残基时，原子名称改动出了bug了？一步的打算是重新上述步骤，在生成rtp文件后插入时，再重新对非标准残基命名规范一下，如果实在解决不了，就只能上传附件求助了

lisanoid

7. 继续运行pdb2gmx ，修复bug
残基SEP的操作步骤同前，操作完以后，再运行6.pdb2gmx命令，此时出现提示，有34个atom 信息missing，经分析后认为，是处理过程中TPO和SEP的H原子未做处理导致的，此时有两个处理办法
一是参考李继存老师博文，他提到的处理方法原文为“（4）如果你需要使用pdb2gmx -ignh来自动添加氢原子, 那就要在相关分子类型的hdb文件中创建新的条目, 指定添加氢原子的个数个方式. 或者, 你也可以新建一个hdb文件, 使用与默认不同的文件名, 以避免直接修改力场自带的hdb文件.”
二是观察PDB源结构中，配体和受体结合情况，如果小分子配体距离比较远，一般分子可能不会与这几个非标准氨基酸产生非键作用（尤其是氢键），可以考虑忽略。
在discovery studio中观察发现，配体主要结合的氨基酸在350-480号残基之间，因此暂不处理氢原子，在命令中加入missing
gmx pdb2gmx -f protein.pdb -o protein.gro -p topol.top -missing
此时gromacs成功运行，并生成protein.grohe topol.top文件。生成了3个肽链的文件。
观察pdb结构发现，这3个肽链结构是重复的，因此后续我会只保留1个肽链，生成拓扑文件后跑MD。观察结果。