计算化学公社
标题: 求助关于Gromacs模拟蛋白质的一些小白问题 [打印本页]
作者Author: 十万嬉皮 时间: 2020-11-1 20:42
标题: 求助关于Gromacs模拟蛋白质的一些小白问题
下面的问题我基本都在论坛里搜过了，但是可能有部分没搜到...如果存在有说明的问题，希望老师能给我个链接。感谢各位老师
1. 第一个问题是，在利用Autodock做完分子对接后，到底什么样的结合模式才应该拿来做分子动力学分析呢。如下图所示，一张是在c端和n端中间夹着的，第二张则是完全被包裹在蛋白内。这个蛋白是一种昆虫的气味结合蛋白（odorant binding protein），含有3个二硫键和6个alpha螺旋的小分子蛋白质；
2. 第二个问题是我在将蛋白和配体拆分的时候，将配体部分截取出来，结果放在gaussion view里面不会自动加氢。这个问题我之前在论坛里面搜如何让gview自动加氢，虽然知道了添加自动加氢的位置在哪里，但是对于这个问题的答案是没有用的。但是就在最近几天，因为配体的问题我跑去acpype的在线网站去做配体的itp文件出现了很多麻烦，于是我一直怀疑是自己在这一步出了问题，之后惊奇的发现，把配体的pdb文件最后一行的H原子都删掉，再放进guassionview里面，竟然全部自动加上氢了，但是很奇怪的是今天跑完rmsd，明明mdp文件里设置了消除体系平动，但是蛋白和配体还是一起动了，我一开始是怀疑自己是不是在设置索引文件的时候出问题了，不小心选错了别的，但是我重新好好设置索引文件后发现并没有出问题。这和以往很不一样，以前都是按照这个流程走的，不会出现平动现象，于是我认为可能会和我删除这个H原子有关联，现在正在重新跑（配体文件上传）；
3. 第三是我想说明acpype网站的机制，如果是小分子在queued的时间特别长，一定要注意是不是自己的配体有问题。因为看教材，在苄脒阳离子示例中，有“-n 1”这个添加项，所以这次我就在net charge中选择了1，然后就一直处于排队状态中，因为以前都非常快，最慢也不会说超过10分钟...而这次排队就一直排了两天，于是我重新提交了新的任务，重新注册账号，提交新任务，即便是改正过哦后的也不行，最后干脆决定自己用ambertools来做了，结果我发现自己加“-n 1”也是报错的，很无奈我四处求助，最后我试了一下去掉了这个选项，竟然成功了...问了sob老师是可以用的。于是当晚我立马就把之前的任务全部取消，提交了一个新的小分子，没有把net cahrge调高，果不其然，提交上去就进入了running模式，很快啊！结果就出来；
4. 第四个问题是昨天我按照教材里给的md文件，只改了步数为1亿步（步长2fs），想跑200ns的md分析，结果出来后还是出现了蛋白的平动现象，之后我把md.gro载人了vmd发现蛋白已经靠近了水盒子的边缘，之后在群里问sob老师，我就把盒子在vmd中使用 pbc box 画了出来，每隔50帧显示一次（应该是50，stride我填的51..），结果发现果然蛋白一部份区域撕裂的很厉害，之后我用通过索引文件用trjconv修复了protein_lig并得到了system的轨迹文件，protein_lig的确是不会出现rmsd的突然上升状况，但是配体却出现了，或者有可能是我的配体没有达到稳定（下图显示的是配体（ligand/MOL）rmsd）；
(, 下载次数 Times of downloads: 48)
5. 最后一个问题是关于rmsd怎才能算是相对平衡或者说是稳定呢..之前也有问过sob老师，但是sob老师说是稳定的（不是上面的图）。但是其实我还想就是自己学一下这个到底怎么才能算是稳定。就比方说相对于backbone的蛋白，对于backbone的配体，和对于complex的complex，之前好像搜到过一个帖子，说是上下波动不超过1埃就算是相对稳定，是这样吗，或者说是有更加贴切具体一些的说法或者证明...
如果各位老师有空可以稍微看一下，我自己知识量太少，所以很有可能犯一些很低级的错误或者有一些问题描述的比较低级，希望各位老师多指正一下。
作者Author: 十万嬉皮 时间: 2020-11-1 20:47
蛋白的图片不知道为什么没有上传上去。
作者Author: sobereva 时间: 2020-11-1 23:31
怎么正确贴图在置顶的新社员必读贴里已经说得很明确了，这里也强调了http://bbs.keinsci.com/thread-18961-1-1.html
作者Author: 十万嬉皮 时间: 2020-11-2 08:57
对不起sob老师，我下次发帖子一定记住
作者Author: sobereva 时间: 2020-11-2 13:04
1 都可以跑动力学，看看结合的动力学的稳定性，算算结合自由能
2 注意正确拼写GaussView
消除体系平动和消除蛋白+配体复合物的平动是两码事。复合物移动了说明消除整体运动的组没定义对，跟H毫无关系
3 阳离子显然要-n 1，中性显然什么也不用写。每个选项都要搞清楚意义何在。诸如我的gromacs培训班上对乙醇用acpype产生拓扑文件就没加这个
4 这种事必须结合轨迹来看到底发生了什么
5 肉眼看上去，RMSD曲线最后阶段没有整体变化趋势就行了。
作者Author: 十万嬉皮 时间: 2020-11-2 20:10
抱歉老师回复晚了！
谢谢sob老师的逐一解答，尤其是第一条，今天跑完了新的动力学分析，蛋白配体复合物的rmsd非常的稳定，蛋白和配体也分别很稳定。
但是自由结合能里面最后一个氨基酸（第118个）的贡献值时60kj/mol，我看很多文献里面关于我这个蛋白的基本上都是中间的两三个氨基酸（例如第30-70内的某几个）贡献的结合能较高（一般都是1~2kcal/mol左右），但是我这个跑完后非常担心是不是在Autodock里面选取的构象不好，因为在Autodock里面在最后的Analyze-conformation-play-(倒数第二个符号）选取Build h-bond，发现这个配体一直和末端的氨基酸产生氢键，于是自己看了关于Autodock更多一些的教学，看到Autodock vina的作者在官网发布的教学，发现gridbox只选取之前文献报道的坐标位置和大小，于是自己就把gridbox调小了，发现可以和中间的一些氨基酸产生氢键了，但是我觉得应该尊重一下事实（可能还是我操作有问题），我决定把末端突变一下去做实际的竞争实验看看到底有没有影响，再跑一下和中间几个氨基酸产生氢键的动力学分析，算一下自由结合能。
还有谢谢sob老师纠错
，一定改正。
作者Author: 十万嬉皮 时间: 2020-11-3 20:58
【关于后续的分享】1. 蛋白建模：由于Modeller建模的问题耽误了很久的时间，因为以前一直是用swiss-model（https://swissmodel.expasy.org/）建模，但是经常会出现缺斤短两的现象（例如我的氨基酸序列比模板的氨基酸序列会少两个，他竟然会自动去掉，这就让人很苦恼，因为后续做各个氨基酸自由结合能的贡献值，如果这两个很关键，那不就出大问题了...）。但是Modeller的单模板建模经常会出现一个问题，在进行SAVE的Verify 3D，ERRAT，拉式构象图的检测中，即便是经过loop循环的优化后，也会打不过swiss-model的建模，到最后我决定在ERRAT氨基酸错误（红色）变成警告（warning）后就拿去用一下；
2. 分子对接（Autodock）：（在这里纠正我之前的一个错误，那个最后一个是 -60 kj /mol 的所谓 “氨基酸” 其实是配体，今天在算完MM/PBSA后发现energy_Mapln文件中排在最后的"MOL"字样，才发现自己是真的粗心大意，错把配体当成氨基酸了。）在这次对接的时候，我刻意的选取了在结合腔内的部分作为Gridbox，然后跑完程序发现一个问题，就是在最后显示结合能的阶段，以往这个蛋白和配体的结合能都是 -4.x，在我人为的干扰之后，就变成了 -1.x，以前师门的师姐说过，这个值越低越好，所以有点担心，但是有一点就是在这个构象当中这个配体不再和蛋白末尾的氨基酸形成氢键，这个让我感觉可以一试，于是我就继续进行下去了；
(, 下载次数 Times of downloads: 65)
3. 分子动力学： 今天我尝试了2ns（基本无意义）和20ns的MD，结果有点让人诧异，20ns的数据，蛋白、配体、蛋白-配体复合物三者竟然在2~3ns之后就趋于稳定了，我出于好奇，继续计算了自由结合能，结果也是出人意料地比较好...让人有点措手不及，但是有一点还是印证了我的想法，就是即便是我把Gridbox人为的调整，使配体和蛋白末端氨基酸不其反应，但是最后还是末端的那个氨基酸赢了（就是还是最末端的氨基酸贡献的结合能最高，虽然不是最后一个，但是是倒数第三个）
(, 下载次数 Times of downloads: 55)
(, 下载次数 Times of downloads: 63)
(, 下载次数 Times of downloads: 56)
(, 下载次数 Times of downloads: 55)
4. 今晚最后做一次时间稍长一些的（200ns）动力学分析，看看结果。
作者Author: 十万嬉皮 时间: 2020-11-4 22:30
1. 今天跑完了200ns的动力学之后，发现RMSD的确在大约100ns左右的时候，有略微的攀升，并且最终在约3.5埃左右的范围内波动，但是如此我便有了一个新的问题，那就是如果我只跑前100ns的话，是否也会出现稳定的情况，这样的话可以得到认可吗？如果我跑在200ns往后的动力学的话，是否会继续出现RMSD的升高？
(, 下载次数 Times of downloads: 55)
2. 之后进行了氨基酸残基的自由结合能计算（取150-200ns -dt 100），发现了其中一个小问题，就是相比较前20ns的情况，有一个氨基酸残基的自由结合能低于 4 kj/mol 了。之后又取了 125-175ns （就不付图了）的轨迹做自由结合能计算，-dt 100 发现结果是不同的，不过这个也可以理解，我认为选取的帧数和自由结合能计算的结果是相挂钩的，也就是说，如果我取稳定时期的帧数越多，那么可靠性应该越高（大概哈哈）
(, 下载次数 Times of downloads: 53)
作者Author: k64_cc 时间: 2020-11-5 17:25
大分子采样还是用Parallel Tempering吧，除非跑超长（us level）的MD，不然brute force MD总是有被怀疑采样不充分的风险。
作者Author: 十万嬉皮 时间: 2020-11-5 23:08
抱歉老师回应晚了...这个蛋白分子量比较小，只有120aa左右的氨基酸，用的mdp文件是根据sob老师教学材料里提供的，因为3ATL又240aa左右，所以我就直接拿去用了，这个还需要换成并行退火吗
作者Author: k64_cc 时间: 2020-11-6 00:24
本帖最后由 k64_cc 于 2020-11-6 00:28 编辑
倒不至于谈需要不需要。现在的情况是你用得起，那就用呗，肯定比不用强啊。
用与不用和分子量关系不是特别大，和你要跨越的能垒高度关系比较大。我之前跑10个碱基的DNA序列，用PT和不用都能看出显著差别来。建议还是用了吧，没啥坏处。
作者Author: 十万嬉皮 时间: 2020-11-6 17:53
谢谢老师指点，那么这个只需要在mdp文件中更改设置吗...还是需要其他的操作，谢谢老师回答
作者Author: k64_cc 时间: 2020-11-7 18:08
我记得官网给了教程，你找找……
作者Author: 十万嬉皮 时间: 2020-11-7 22:16
好的，谢谢老师，我去看下
作者Author: 谭谭野野 时间: 2022-1-11 10:03
能加个微信吗？有问题向您请教，微信号t519136647
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |