请问MD优化docking结果时，如何从MD轨迹得到pose呢？gmx cluster似乎不适用于这一场景

papercup

请教各位老师们大佬们，

我目前在疑惑如何用MD做docking后小分子pose的refinement才是合理的，这主要又分为以下两个问题：
（1）怎么从production MD轨迹得到最终refinement后的小分子pose？
（2）如果（1）的答案是聚类，gmx cluster能做这件事吗？gmx cluster以同一个group消除平移旋转和算RMSD，似乎不适用于这一场景


===== 问题详述 =====


（1）production MD的输出是一段轨迹，但refinement的结果应当是一个单一的pose，那如何从一段轨迹最终得到一个改进后的小分子pose才比较合理呢？

用MD去做docking结果的refinement，应该也是一种普遍的做法了，因此我也查阅了一些文章，但是感觉这些做法中，得到最终小分子pose的方式都不太靠谱。
https://pubs.acs.org/doi/10.1021/acs.jcim.0c00057
像这篇文章，（虽然也拿轨迹的RMSD排除了一些decoy）读下来大意是取了轨迹最后一帧作为最后输出的pose，感觉不太靠谱。
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6637856/
这篇文章的话，大意是先MM能量最小化，再MD一段时间，取MD最后一帧再MM能量最小化，得到最后的pose，感觉也不太靠谱。

这个问题其实也不限于小分子pose，用MD优化蛋白质结构应该也需要解决从轨迹到结构的问题吧，不知大家都是如何做的呢？

询问同事并论坛搜索后，得到的答案是做聚类，这个感觉就靠谱多了。
http://bbs.keinsci.com/thread-7566-1-1.html
但此时就遇到了问题（2）。



（2）gmx cluster聚类是以同一个group消除平移旋转和算RMSD，是否意味着不适合提取小分子的代表性pose呢？

运行gmx cluster时会有类似这样的输出：

1. ......
   
2. 
   
3. Command line:
   
4.   gmx cluster -f ../md.xtc -s ../md.tpr -cutoff 0.25 -b 500 -wcl 10 -method gromos -sz -n ../index_protein_ligand.ndx
   
5. 
   
6. Using gromos method for clustering
   
7. Reading file ../md.tpr, VERSION 2019.1 (single precision)
   
8. Reading file ../md.tpr, VERSION 2019.1 (single precision)
   
9. 
   
10. Select group for least squares fit and RMSD calculation:
    
11. Group     0 (         System) has 46319 elements
    
12. Group     1 (        Protein) has  4405 elements
    
13. Group     2 (      Protein-H) has  2214 elements
    
14. Group     3 (        C-alpha) has   272 elements
    
15. Group     4 (       Backbone) has   816 elements
    
16. Group     5 (      MainChain) has  1090 elements
    
17. Group     6 (   MainChain+Cb) has  1341 elements
    
18. Group     7 (    MainChain+H) has  1352 elements
    
19. Group     8 (      SideChain) has  3053 elements
    
20. Group     9 (    SideChain-H) has  1124 elements
    
21. Group    10 (    Prot-Masses) has  4405 elements
    
22. Group    11 (    non-Protein) has 41914 elements
    
23. Group    12 (          Other) has    27 elements
    
24. Group    13 (            MOL) has    27 elements
    
25. Group    14 (             CL) has     4 elements
    
26. Group    15 (          Water) has 41883 elements
    
27. Group    16 (            SOL) has 41883 elements
    
28. Group    17 (      non-Water) has  4436 elements
    
29. Group    18 (            Ion) has     4 elements
    
30. Group    19 (            MOL) has    27 elements
    
31. Group    20 (             CL) has     4 elements
    
32. Group    21 ( Water_and_ions) has 41887 elements
    
33. Group    22 (    Protein_MOL) has  4432 elements
    
34. Select a group: 22
    
35. Selected 22: 'Protein_MOL'
    
36. 
    
37. Select group for output:
    
38. Group     0 (         System) has 46319 elements
    
39. Group     1 (        Protein) has  4405 elements
    
40. 
    
41. ......

复制代码

gmx cluster用来least squares fit和RMSD calculation的group是同一个。
但在分子对接refinement的应用场景下，此处应当以Protein或Backbone消除平移旋转，再以MOL计算RMSD才比较合理吧。
因为优化docking结果的小分子pose时，显然是不能去对齐小分子也不能忽略小分子的平移旋转的。
若建立Protein_MOL的group，则RMSD会受Protein大小的影响，使得cutoff很难确定，带来另一方面的麻烦。

不知道gmx cluster支持这样的用法（以蛋白质对齐，再以小分子计算RMSD）吗？
还是说只能导出轨迹后再手动实现聚类呢？

panernie

浅谈PCA与g_covar+g_anaeig+ddtdp+sigmaplot做自由能面图的方法-http://sobereva.com/73

sobereva

“以蛋白质对齐，再以小分子计算RMSD”是可以实现的
gmx cluster命令可以加上-nofit，这样算RMSD前就不会先做fit。你可以先用trjconv依据蛋白质进行fit，然后gmx cluster的时候只算小分子的RMSD

papercup

sobereva 发表于 2021-9-1 21:24
“以蛋白质对齐，再以小分子计算RMSD”是可以实现的
gmx cluster命令可以加上-nofit，这样算RMSD前就不会 ...

非常感谢sob老师的解答

今天按照您回复中的做法尝试了一下，果然成功了，十分便捷。
这下就省事多了，先前我还想着可能需要用mdtraj分析轨迹+自行实现聚类才能做到，想不到其实只需合理使用gmx cluster就行了。

十分感谢！

papercup

panernie 发表于 2021-9-1 21:14
浅谈PCA与g_covar+g_anaeig+ddtdp+sigmaplot做自由能面图的方法-http://sobereva.com/73

感谢大佬回复

文章中的做法，我理解是PCA降维至二维后，再人工地去挑选最大的cluster中的帧
但其实我这个问题里，各pose之间的RMSD能很快速且精确地求出来，直接用RMSD为距离做聚类就行

今天我尝试了sob老师回复中的做法，直接用gmx cluster以RMSD聚类就能很快出解了，能免去降维和人工去看的麻烦

frank666

您好，使用MD来得到对接分子refinement的过程是我比较感兴趣的研究方向，看了您的帖子觉得很受启发，但也产生了一些遗憾。
1.关于您说“以蛋白质对齐，再以小分子计算RMSD”，在gmx cluster中不是可以先选Protein再选MOL吗，得不到一样的效果吗？
2.gmx cluster加-nofit命令与不加-nofit我都试了一下，两个命令针对相同体系，我都是重复选择两次MOL，-nofit命令下得到的RMSD范围要大的多，请问fit的本质，到底是什么呢？为什么得到RMSD范围差那么多呢
3.因为我也想根据您的思路去得到分子对接结果的refinemen，根据您“以蛋白质对齐，再以小分子计算RMSD”意思，请问是先在trjconv -fit rot+trans，选择两次蛋白组，然后在gmx cluster -nofit，选择两次配体组，最后根据cluster的结果选取对接构象吗
无论您是否回复，都感谢您提供了如此精彩的帖子，您研究思路的严谨性、科学性都让我深感佩服。因为我还入门不久，如果以上问题让您倍感无趣，请您原谅，如果您能抽空回答，我将感激不尽。