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

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后续及结果分析
   本帖拟分步记录处理过程，由于本人水平较低，处理过程出现了问题，正好请社长和多多指教。

lisanoid

   这次研究的目标蛋白是PRKC  蛋白质pdbID：3IW4。从RCSB网站（https://www.rcsb.org/）下载结构，
  采用Discovery Studio裁剪掉水分子，将配体剪切后新生成ligand.mol2文件，剩下的蛋白质结构保存成protein.pdb文件，此时如果对该pdb文件做pdb2gmx，出现报错fatal error，提示 residue TPO638 is of type 'Other'.....
   用文本编辑器打开pdb文件，发现有missing residue和missing atoms，为防止报错是由于missing残基导致，将蛋白质的pdb文件在pulchra中处理，处理后的文件运行pdb2gmx后，仍然出现报错，目标明确为638号残基TPO。此时也能生成拓扑文件topol.top，将该文件保留，后续打算都跑一下MD做个对比。
    明确任务后，先做第一步裁剪残基这一步，用disovery studio 这个软件操作简单傻瓜化，还有VMD应该也能实现，或者直接用文本编辑器应该也可以，但我尝试操作过但没有成功。剪切得到的TPO文件再用gaussview打开，在N端接乙酰基，在C端接甲胺基，检查分子中原子价键，连接方式以及基本的结构没有错误后保存。
     本步主要任务是明确错误目标，其次是构建结构正确的残基做计算。通常来说，做生物的很多人化学基础薄弱，对化学结构并不敏感，这一步错了后面处理过程一点意义也没有，因此需要仔细检查。
      得到结构以后在gaussview中保存为tpo.gjf文件，用文本编辑器把gjf文件打开，把第二行# hg3.....替换为以下关键词： # opt freq b3lyp/6-31++g geom=connectivity  保存gjf文件后，用gaussian优化结构。（这一步我在做的过程中出错很多，在参考社长文献时，使用# B3LYP/6-311G** em=GD3BJ opt scrf=solvent=water作为关键词，报错l19999，说是有空格问题，一直没搞清楚。这步操作需要有一些量子化学基础。）
     普通的4核或者8核笔记本电脑，这一步操作比较耗时，后续操作我们周一继续。
      
   

lisanoid

2. 计算RESP电荷
  上述结构优化完以后，得到TPO.out文件。我的操作是先用win版的gaussian，然后utilities选项下生成formchk文件，选择TPO.out文件，高斯就生成了TPO.fch文件，此文件带有multiwfn能识别的波函数数据。根据http://sobereva.com/441  还有http://sobereva.com/637等等介绍可以先学一下。
    计算RESP电荷，具体做法是，将生成的fch文件导入multiwfn，然后进入Multiwfn主功能7的子功能18 ，再用选项1 计算标准RESP电荷，导出chg文件后，最后1列就是RESP电荷。
   本步和第1步可以联合解决，具体做法可以用http://sobereva.com/637 的方法，直接调用懒人脚本。如果为了发文章而又没有授权的量子化学软件，建议参考这一步，调用orca解决。这样的话，需要安装虚拟机系统，安装orca和multiwfn的linux版本，社长博文中都有详细介绍，可以自行学习。
     如果再懒一些，我觉得这一步也能省略，因为常规的药物分子对接和虚拟筛选，使用最广泛的免费软件autodock，默认的原子电荷是Gasteiger，精度上与resp根本没法比。所以如果这一步不更新电荷，后续操作做出来的MD应该也比分子对接的结果要好。

lisanoid

3.构建rtp文件
    根据sobtop的FAQ5 ，rtp文件生成方法为：然后在Sobtop里载入模型分子的mol2文件，从chg文件里载入原子电荷，然后选产生rtp文件，其间让Sobtop自动指认AMBER原子类型，并要求从预置力场库文件里取合适的参数，缺的让Sobtop基于Hessian算出来。然后把得到的rtp里的字段恰当处理从分子状态变成为残基状态后，挪到GROMACS力场目录下的rtp里。
    以下是操作步骤：
将优化结构的.out文件用gaussview打开，生成mol2文件。用sobtop载入文件，
选项7 assign 原子电荷
选项10采用4计算的chg文件
把第2步生成的chg文件拖入程序，此时resp电荷被加载
0 return
3 产生rtp文件
4指认amber原子类型
0返回 //此时已经加载参数
4 Same as option 3, but missing ones are arbitrarily guessed
回车 文件生成在sobtop文件夹下
0 退出
此时得到rtp文件

lisanoid

4 rtp文件的编辑
rtp文件生成后，不能直接使用，因为当时生成时是加了甲胺基和乙酰基的，而正常的肽链中，端基氨基酸少甲胺基或乙酰基，而非端基氨基酸两个都没有，因此我们需要将rtp文件中涉及甲胺基和乙酰基的统统删掉，同时我们在做rtp文件时，原子编号是重新生成的，与pdb文件中的原子名是不一样的，因此我们还要重新编辑rtp文件，使之能够在放入gromacs以后能被正确调用，生成所需的蛋白质拓扑文件。本步操作繁琐，需要耐心对照化合物结构一点一点对应完成。
此步骤的操作方法如下： 将第3步中生成的TOP.mol2分子文件用gaussview打开，在view选项中选择label，观察甲胺基和乙酰基中各个原子的编号，我生成的文件中是从N20到O31，把20号原子之后的在rtp文件中统统删掉，另存rtp。
再打开蛋白质pdb文件，将pdb文件中的原子名（第三列，原子序号之后的那一列），按照一一对应的原则跟之前rtp文件中的原子名改成一样的。

lisanoid

该步操作忽视了sobtop产生的rtp是不支持gromos力场的，因此拷入gromos54a7是错误的，此步可以拷贝到amber相应力场的ff文档中，如amber99sb-ildn。
5.把rtp文件整合进gromacs中
得到的TPO.rtp文件，需要在gromacs运行时被正确的调用，因此需要把相关的信息都写入gromacs相应文件中，李继存老师博文https://blog.sciencenet.cn/blog-548663-1075988.html中，对此有详细描述，内容比较多，参考内容一步一步完成。
我的操作如下：
--把gromacs文件中top\gromos54a7.ff文件夹拷贝一份保存好，此步骤繁琐易错，为防止出问题备份好原始文件；
--把tpo.rtp文件中的信息，按照sobtop程序给出的提示分别转移，将生成的TPO.rtp文件放到54a7.ff文档中，
--按要求将原子类型信息挪到atomtypes.atp文件夹中，在rtp文件中atomtypes部分，拷贝到54a7.ff文档中ffnonbonded中，整理格式。
--把生成的rtp文件中，[MOL]名字头改为[ TPO]，rtp文件拷贝到54a7.ff文件夹下，并将rtp文件中开头的两节信息：atomtypes和原子类型（？此处在rtp文件中生成后无标题，我也不知道是啥）删除掉，此时rtp文件内容第一部分应为bondedtyoes。
--在../top文件夹中的residuetypes文件中，加入[TPO],定义为protein
把上述更改的文件整合进top\gromos54a7.ff文件夹中。

此时可以调用该残基信息了。

lisanoid

6.运行pdb2gmx 修复bug
把protein.pdb文件拷贝到合适位置，运行程序
gmx pdb2gmx -f protein.pdb -o protein.gro -p topol.top
此时出现新提示，657SEP也是other，意思就是这个也需要单独处理。太可恶了，上次报错提示了个638号氨基酸 ，不一次性全提示完，这样又要重复上述步骤再建立一个新非标准氨基酸。
这次有了经验，继续按原方法优化氨基酸结构：
    剪切出残基SEP，加甲胺基，加乙酰基，然后用懒人脚本，调用multiwfn+orca，优化结构，计算RESP电荷。
同时检查了后续氨基酸，确定再没有非标准氨基酸了，等结构优化以后再重复上述操作。

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。观察结果。

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

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

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

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

lisanoid

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

hongyi166_

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

liu840397562

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

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