使用xtb进一步优化分子对接结果，表现上却还不如经验公式，做法是否有误？

papercup

请教各位老师们大佬们，我目前在寻求一种方法对蛋白质小分子对接（protein-ligand docking）的结果做进一步的优化（refinement）。
搜索论坛后，参考了一些帖子的做法，比如这个：
http://bbs.keinsci.com/thread-19956-2-1.html

马上用 xtb 试了试，但批量跑了几百个 case 后，发现在 docking 的表现上还不如一些经验公式（比如 autodock vina 自带的能量打分函数）好。
请问是否是我对 xtb 的使用方法出错了？

===== 分割线 ======

我目前的大致做法是：
1、取几百个结晶的实验结构，直接以实验结构为起点开始优化，如果 xtb 比较准，那优化结果应当与起点偏离地很小；
2、做截断，只取距离 ligand 5A内的所有 protein residues ，并补氢；
3、将 residues 和 ligand 组成 complex，总体上约有700～1000个原子；
4、fix 住所有 protein 部分，用 xtb 默认设置（GFN2-xTB）对 complex 进行能量最小化；
5、计算 ligand 优化结果与起点的RMSD，若RMSD>2则认为偏离了起点。

批量跑了几百个 case 后，发现有约 10% 的 case 都偏离了起点。
但若以 autodock vina 自带的能量打分函数优化实验结构，100% 的 case 都可以停留在起点。
类似地，我也尝试了批量用 xtb 进一步优化 vina 输出的 docking 结果，反而有部分原本正确的 docking 结果被 xtb 优化到了错误的位置。整体上，正确的 case 数也更少了。


我的疑问是：
（a）xtb 作为半经验的量化计算程序，理应在能量计算上比传统力场和经验公式更加准确，为什么在 protein-ligand docking 这个问题上，表现甚至还不如传统的 autodock vina 经验公式呢？
（b）我的 xtb 使用方式是否存在错误？
（c）有没有什么已验证可行的 refinement 方式能推荐一下呢？基于力场/量子化学的都可以。

wzkchem5

优化时考虑水了吗？残基的质子化情况都验证过吗？有没有正确考虑蛋白整体的电荷，或者加抗衡离子中和蛋白的电荷？

sobereva

GFN-xTB优化精度本来就比较有限，单分子尚且如此，更不用说相互作用繁多、势能面复杂的蛋白质口袋区域了。虽然打分函数经验性更强，但那毕竟是为蛋白质-配体结合问题特意设计的，所以结果更好也并不意外。

你的计算用的模型的合理性也对结果有明显影响，细节没有在你的描述中体现。对于发现有问题的那些体系，最好再用更靠谱的量化级别优化，比如ORCA里的B97-3c，看看是否能得到你期望的结果，弄清楚到底是GFN-xTB的问题还是模型构造的问题。

不要把蛋白所有原子都冻住，最好让侧链可以活动，这样的话和配体相互作用的残基会自发弛豫，增强相互作用，可能令配体原本待不住的情况最终能近似呆在原位。另外，本身X光衍射的精度也是有限的，允许适度弛豫会更合理。另外你没有说你取的实验结构的解析度是多少，只有解析度较高的才能用于对比测试。

fhh2626

蛋白质-配体相互作用熵贡献很明显，打分函数在拟合过程中应该是考虑了一部分熵贡献的，但是xtb算单点或者优化是不包括这个贡献的

papercup

wzkchem5 发表于 2021-6-30 00:15
优化时考虑水了吗？残基的质子化情况都验证过吗？有没有正确考虑蛋白整体的电荷，或者加抗衡离子中和蛋白的 ...

非常感谢大佬的回复

1、优化时没有显式地加水，只是运行xtb时加上了 --gbsa

2、
质子化我的做法其实不太好，我做截断时，只是基于修复了侧链和non-standard residues后的 protein.pdb 文件简单地选择了距离较近的残基。删去所有氢的行后再用 openbabel 补了氢
这样做可能切断了 -C-C- 以外的键，也切断了一些环，而且原先是 -OH 的地方现在可能被错误地补上了 -H
如果只是氢的个数的话应该是正确的，有可视化地看过，从 xtb 输出来看体系的电荷也都是 0.0

3、
我今天看了一下，截断后的体系整体的电荷，从 xtb 的输出文件上看，均为 0.0
但如果运行 gmx pdb2gmx 去查看最初的输入的 .pdb 蛋白质文件的电荷，大部分蛋白质都是带一定的正负电，不为 0.0 的
不知这是否会带来很大的影响

4、因为没有显式加水的关系，好像也无法通过替换溶剂添加离子中和电荷

papercup

sobereva 发表于 2021-6-30 03:00
GFN-xTB优化精度本来就比较有限，单分子尚且如此，更不用说相互作用繁多、势能面复杂的蛋白质口袋区域了。 ...

非常感谢sob老师的回复

我其实最担心的也是我的计算用模型的合理性，因为我在这方面的经验很少，完全不能信赖
就像我上一条回帖写的，目前做截断的方式可能也存在一些问题，切断了 -C-C- 以外的化学键，也切断了环。
我现在的想法是更多地去借助一些 MD workflow 相关的第三方工具（比如说 biobb）去完成体系的准备，避免自己去调用工具修复侧链或是修复non-standard residues，一来这样肯定比我个人的 workflow 靠谱；二来我这个 docking 问题本身也需要更自动化更批量的打分方式，若需要大量的手工操作就意义不大了

用准备 MD 体系的方式（但也做截断）准备计算用模型，以 Gromacs 验证这个准备的方式正确后，再将去掉显式水的体系通过 xtb 和您提到的 ORCA 里的 B97-3c 优化。
这样的思路是否是合理的呢？

另外关于解析度，我用的是 pdbbind 数据库下的一部分 case，解析度基本能达到 <2A

papercup

fhh2626 发表于 2021-6-30 11:14
蛋白质-配体相互作用熵贡献很明显，打分函数在拟合过程中应该是考虑了一部分熵贡献的，但是xtb算单点或者优 ...

非常感谢大佬的回复

请问如果是同一个 ligand 不同 pose 的情况，熵是否也会带来影响呢？
因为 autodock vina 的能量打分函数表现上很不错，但从每一项来看好像也并不包含熵，vina的打分是包含5项：
1、范德华力
2、repulsion防止原子过近
3、氢键
4、疏水性
5、一个扣减，若 ligand 内部的旋转轴越多，则评分越低（有一点类似熵的概念，但对于同一个 ligand 来说，这一项是个定值）
在挑选 pose 的时候，vina似乎不考虑熵，只用到1～4的能量项也达成了很好的效果

sobereva

papercup 发表于 2021-7-1 00:18
非常感谢sob老师的回复

我其实最担心的也是我的计算用模型的合理性，因为我在这方面的经验很少，完全 ...

合理

wzkchem5

papercup 发表于 2021-6-30 17:17
非常感谢大佬的回复

1、优化时没有显式地加水，只是运行xtb时加上了 --gbsa

我觉得一方面显式加水可能比较好，不过--gbsa可能也够了，毕竟MMGBSA也属于靠谱性至少不亚于一般docking打分函数的方法，说明GBSA应该不算太差。
另外，电荷的正确性还是需要仔细确定。你看到xtb输出电荷为0，那是因为当命令行不指定电荷时，xtb默认电荷为0，xtb不会自动为分子的各个原子指定形式电荷，这也是所有半经验方法的共性。所以必须你人工确定（截断过的）蛋白，以及（截断过的）蛋白-配体复合物的总电荷分别是多少（只需要知道总电荷就够了），然后用命令行参数告诉xtb。如果实际蛋白电荷不为0，但是你没有指定电荷，那么结果会定性错误。

fhh2626

papercup 发表于 2021-7-1 00:19
非常感谢大佬的回复

请问如果是同一个 ligand 不同 pose 的情况，熵是否也会带来影响呢？

打分函数中的项虽然说是能反应什么作用，但是其实不具有明显的物理意义。vina（或者更新的vinardo）的打分函数是拟合结合自由能得到的，所以熵贡献已经在拟合的过程中隐含的包含在里面了。对同样的ligand不同pose的对比，熵贡献应该小一些，但是也不能肯定就没影响了

papercup

wzkchem5 发表于 2021-7-1 03:42
我觉得一方面显式加水可能比较好，不过--gbsa可能也够了，毕竟MMGBSA也属于靠谱性至少不亚于一般docking ...

非常感谢大佬指出这个错误。xtb需要手动指定总电荷我却都给设成0了，那我之前的测试可真是错的蛮离谱了

昨天今天我都在调试算蛋白质电荷的程序，一种方式是直接数氨基酸 LYS ARG HIS ASP GLU 来算；一种是把截断的蛋白质补好 -OH 和 -H 后用第三方程序去算。要是我补 -OH 和 -H 补得对，电荷算得也对，那两种方式的结果应该会相同。现在也只剩一小部分 case 程序跑完电荷对不上了。
电荷的计算问题解决后，我马上重跑 xtb 测试一下结果。