利用AutoDock进行水分子参与的分子对接

水分子在配体-蛋白质相互作用中的作用至关重要，AutoDock开发小组的研究人员设计了一种新的力场，可以优化水分子参与的分子对接，使得AutoDock程序能够自动预测水分子介导的配体蛋白结合方式，该力场考虑了离散水分子对配体结合的熵和焓贡献，提高了能量估算精度和对接性能。该方法无需事先知道蛋白与水分子的水合状态，在基于结构的药物发现过程中具有较大的应用价值。本文对AutoDock的水合对接方法的精度进行了总结，并且提供了如何利用最新版AutoDock进行水分子参与的分子对接教程。

1 背景

在生理环境中，蛋白质和其他生物结构被水分子包围。当小分子与蛋白质结合时，它必须置换占据结合口袋的大部分水分子。然而，很少有所有的水分子都发生移位。有研究表明，在大多数复合物中(>84%)有一个或多个水存在并与配体相互作用，介导其与蛋白质的相互作用。

有些水分子与蛋白质之间存在强相互作用，或在相似的蛋白中是保守的、它们被认为是靶标蛋白结构的一部分，起着改变结合位点的形状的作用。经典的例子有HIV-1蛋白酶PR和乙酰胆碱受体，其中稳定水的作用是增加抑制剂的亲和力或促进药效团的定义。稳定水可以通过向溶剂中释放有序水而产生的熵增益以提高配体的亲和力。这些策略已成功应用于设计小柱孢酮脱水酶抑制剂和PR的循环尿素抑制剂。弱结合水更容易根据结合配体的性质发挥不同的作用。事实上，同一种水可以被一种配体稳定而被另一种配体置换，如同PARP 抑制剂一样。

然而，水置换并不总是导致亲和力的提高。选择哪种水置换对药物设计可以很重要。常用的对接软件包括AutoDock，通常的流程是在对接前先从蛋白质结构中除去所有的显性水，然后以连续的去溶剂化势能的形式使用隐式溶剂模型。对于已知存在一个或多个水的情况，可以通过包含选定的显式水来模拟多种形态的靶标。

可是，对于所有的配体，保持它们在相同的取向会产生偏差，这将是配体与不同水模式结合面临的问题。例如，PR中结构水的存在不允许循环尿素抑制剂正确对接。此外，比较不同水分子替换情况下得到的结果是一项非常重要的任务，因为很难计算水分子替换在蛋白质中的位移所产生的熵贡献。另一个问题是可以存在于晶体结构中的水的分类。通常定义良好的结构水包含在分辨结构中，而弱束缚水更容易被忽略。

随着计算对接作为药物设计中常用的配体筛选方法的兴起，一个快速、准确的多靶点结合位点水化模型的建立对结合能估算和结果准确性至关重要。目前，这个领域已经发展了不同的定性和定量方法来预测蛋白质结构中水分子的存在或位移所隐含的能量贡献。它们可分为基于结构的方法，它们依靠实验确定的水位置(例如来自高分辨率晶体结构)，又或者是水分子的位置由计算方法确定。

基于实验结构确定水分子的位置有几个缺点：首先，晶体结构与配体结合状态会影响与水分子的结合态，水分子可能通过配体的相互作用而被稳定下来；其次，晶体结构的质量受到多方面的影响，包括解析精度、温度因子、不同的解析方法（如NMR无法解析水分子）。

图1. 水合态示意图

在分子对接中，几个程序可以处理配体对接过程中结合位点的显式水，如GOLD、GLIDE、SLIDE和DOCK这类全原子表征方法或FlexX和FITTED这类球形粒子表征方法的对接程序。这些对接程序大部分依赖于预定的和定向的水的位置，但可建模的显式水的数目受到随之而来的自由度数目增加的限制，典型的如GOLD采用的可建模水分子为3个。此外，与水位移有关的能量常被忽略在评分函数，如FlexX 、GLIDE和DOCK程序。FlexX并不要求靶标结构中存在水，但在预测其位置和改善对接结果方面的准确性据报道很小。DOCK和GOLD则采用包含水分子的对接方式，因此性能相对提高。然而，后者的作者推测，包含可转换的晶体学水可以通过减少搜索空间来完善结果，这可能会使与不同配体的de novo对接比实验结构中存在的更不客观的结果。

图2. AutoDock预测水分子在配体周围分布情况

AutoDock小组将原来的力场修改为在对接过程中包含显式置换水，以在不增加过多计算复杂度的情况下提高对接精度和能量估计。研究人员不是将水放置在靶标表面的特定位置，而是在对接前将水附着在配体上，在搜索过程中不断评估它们可能存在的方式。如果水分子的介导稳定了配体和受体的相互作用，它就会保持；否则，就是被剔除。新的方法分别计算了每个模型水对熵和焓的影响。通过对焓贡献和熵贡献之间的平衡进行评估，来区分易于流失的水或者保守水分子。

图3. AutoDock计算水分子替换引起的能量变化示意图

为使得本方法适用于药物设计和虚拟筛选中的潜在用途，需要保证力场必须是普遍适用的，并且不依赖于以往任何关于水存在的知识。因此，它必须复盖广泛的配体尺寸(从小碎片到药物相像性分子或肽类分子)和具有不同亲疏水性性质的靶标，涉及任意数量的水(包括没有)。本方法首先通过结合使用一个配体水化模型和网格图来描述靶标与水的亲和力，然后经过测试发现本方法在计算资源的需求上也较小，可以满足大批量的筛选任务。

2 水合对接精度评估

下面介绍AutoDock开发的新力场参数在197个配体-蛋白质结晶学结构的训练集上进行的测试结果。这些晶体结构其分辨率为2.5埃或更高，包括50多个不同的蛋白质家族。然后在220个配合物和交叉对接实验的测试集上验证了力场。保守水原子的相互作用能也被用来分类为强束缚或弱束缚。对接结果分析和预测已知晶体学水的位置和稳定性，针对4个案例进行了较为详细的讨论。

如图4 所示，在训练集中，标准AutoDock力场对接成功124例(62.9%)，而新水合对接成功143个案例(73.6%)，性能提高17% (p值，0.007)。有趣的是，在123个复合物中，新方法水合对接计算得到的rmsd值均低于用标准协议得到的rmsd值(图4a )。精度的提高在试验集中也得到了体现，其中标准对接和水合对接成功的分别有135例(61.1%)和161例(72.9%)，表现出19%的性能提升(p值0.001 ,图4b )。

图4. 用于训练和测试集的标准对接和水合对接RMSD分布比较图

3 案例研究

结构1uw6包含1个结合态水分子，其中吡咯烷带电基团与Trp147相互作用，阳离子-π和氢键相互作。其中水分子介导Met114和Leu102骨架与吡啶环之间的氢键(图5c)。标准对接和新的水合对接均能以满意的精度再现实验结合模式( 分别为1.18和 0.42埃,图5 )，尽管前者在无水条件下得到的结果显示吡啶环的翻转取向导致氮与Tyr192羟基之间的弱相互作用。水合对接再现了吡啶环的实验构象，准确预测了结晶水的位置。化合物的N-甲基-2 -苯基吡咯啉对接到和合适的位置，显示了新力场的对接性能的可靠性。

图5. 晶体结构1UW6的标准分子对接（a）水合对接（b）比较

4 实战演练教程

AutoDock开发小组将新的适合水合对接的力场整合进了AutoDock4和现在的新版本的AutoDockVina。在试验中，这种方法在预测小片段分子的结合构象方面表现出了改进，如用于基于片段的药物发现。计算流程可以通过这些步骤加以总结：

1）配体上装有一系列水分子(用虚拟原子表示)，这些水分子可能或不可能参与分子间的相互作用。2）然后在对接过程中使用修改后的AutoGrid分布图，在水放置良好时给予有利的评分，如果水与受体重叠则省略。3）最后，对接结果进行分析，只使用保留的水分子重新评分构型。

在本教程中，我们将在乙酰胆碱结合蛋白( AChBP )结构( PDB编号1uw6)中，将一个破碎块度配体(烟碱)与显式水分子对接。上文强调过缺少水分子会对对接性能产生重大影响，导致打分不准确和/或姿态不正确。这里我们将再先上文中的案例。

完成本案例依赖包：

• ADFR
• Meeko
• Vina1.2.x
• RDKit

ADFR安装包地址：

https://ccsb.scripps.edu/adfr/downloads/ 

Vina获取地址：

https://github.com/ccsb-scripps/AutoDock-Vina/releases

Meeko获取地址：

https://github.com/forlilab/Meeko

本教程所用脚本获取地址：

https://github.com/ccsb-scripps/AutoDock-Vina/tree/develop/example/autodock_scriptshttps://github.com/ccsb-scripps/AutoDock-Vina/tree/develop/example/hydrated_docking

4.1 处理受体

图6. 编号为1uw6的蛋白质晶体结构

根据编号获取结果，然后保留A和B链，删除其他的内容，使用Chimera加氢保存为1uw6_receptorH.pdb。使用ADFR的prepare_receptor命令生成受体pdbqt文件

prepare_receptor -r 1uw6_receptorH.pdb -o 1uw6_receptor.pdbqt

4.2 准备配体

图7. 本案例中小分子结构

将配体的三维结构保存为sdf格式文件，使用mk_prepare_ligand.py对其进行加氢以及质子化处理。同时加上-w参数将显示水分子添加进配体文件中。

mk_prepare_ligand.py -i 1uw6_ligand.sdf -o 1uw6_ligand.pdbqt --add_hydrogen --pH 7.4 -w

4.3 生成亲和力分布图

这里需要根据配体的氢键供体和受体原子类型来定义有水分子参与的配体的亲和力分布图。选项-y指定我们要自动地将网格中心设置在配体周围。-p是指定配体原子类型需要包含在GPF文件中，同时忽略W（虚拟原子）。其中OA（氢键受体）HD（氢键供体），手动添加OA原子类型后，应该有一个名为1uw6 _ receptor . gpf的GPF文件，命令如下：

pythonsh <script_directory>/prepare_gpf.py -l 1uw6_ligand.pdbqt -r1uw6_receptor.pdbqt -y -p ligand_types='A,NA,C,HD,N,OA'

GPF参数文件内容：

npts 40 40 40                        # num.grid points in xyz
gridfld 1uw6_receptor.maps.fld # grid_data_file
spacing 0.375                        # spacing(A)
receptor_types A C NA OA N SA HD # receptor atom types
ligand_types A NA C HD N OA # ligand atom types
receptor 1uw6_receptor.pdbqt # macromolecule
gridcenter 83.640 69.684 -10.124     # xyz-coordinates or auto
smooth 0.5                           # store minimum energy w/in rad(A)
map 1uw6_receptor.A.map              # atom-specific affinity map
map 1uw6_receptor.NA.map             # atom-specific affinity map
map 1uw6_receptor.C.map              # atom-specific affinity map
map 1uw6_receptor.HD.map             # atom-specific affinity map * ADD OA IF NOT PRESENT *
map 1uw6_receptor.N.map              # atom-specific affinity map
map 1uw6_receptor.OA.map             # atom-specific affinity map * ADD OA IF NOT PRESENT *
elecmap 1uw6_receptor.e.map          # electrostatic potential map
dsolvmap 1uw6_receptor.d.map              # desolvation potential map
dielectric -0.1465                   # <0, AD4 distance-dep.diel;>0, constant

然后使用Autogrid生成有水分子参与的亲和力分布图

autogrid4 -p 1uw6_receptor.gpf -l 1uw6_receptor.glg

python <script_directory>/mapwater.py -r 1uw6_receptor.pdbqt -s 1uw6_receptor.W.map

检查生成的map文件

ADD PWD AND FILE SUMMARY
  receptor : 1uw6_receptor.pdbqt
      OA map -> 1uw6_receptor.OA.map
      HD map -> 1uw6_receptor.HD.map
 => Water map weight : DEFAULT [ 0.60 ]

  MapWater generator
 =====================
  mode      : BEST
  weight : 0.6
  HD_weight : 1.0
  OA_weight : 1.0
  entropy : -0.2

     Output info
  --------------------
  filename : 1uw6_receptor.W.map
  OA points : 91.73%
  HD points : 8.27%

  lowest map value : -0.99
  highest map value : -0.01

4.4 运行Vina进行对接

上述准备工作完成后，就可以使用最新版的Vina进行水合对接啦

vina --ligand 1uw6_ligand.pdbqt --maps 1uw6_receptor --scoring ad4 --exhaustiveness 32 --out 1uw6_ligand_ad4_out.pdbqt

对接完毕，预测结合自由能为-8kcal/mol

Scoring function : ad4
Ligand: 1uw6_ligand.pdbqt
Exhaustiveness: 32
CPU: 0
Verbosity: 1

Reading AD4.2 maps ... done.
Performing docking (random seed: -655217817) ...
0% 10 20 30 40 50 60 70 80 90 100%
|----|----|----|----|----|----|----|----|----|----|
***************************************************

mode | affinity | dist from best mode
     | (kcal/mol) | rmsd l.b.| rmsd u.b.
-----+------------+----------+----------
   1 -8.077 0 0
   2 -7.63 2.038 2.684
   3 -7.382 2.378 2.747
   4 -7.27 2.063 2.538
   5 -7.138 1.861 5.391
   6 -7.129 2 2.542
   7 -7.078 3.307 5.442
   8 -7.065 2.22 4.872
   9 -7.051 3.135 5.636

使用Pymol打开1uw6_ligand_ad4_out.pdbqt以及受体文件，可看到结合模式如下，与文献报道结合模式基本一致。

图8. 结合模式

参考文献

1. Forli S, Olson A J. A force field with discrete displaceable waters and desolvation entropy for hydrated ligand docking[J]. Journal of medicinal chemistry, 2012, 55(2): 623-638.

2.Forli S, Huey R, Pique M E, et al. Computational protein–ligand docking and virtual drug screening with the AutoDock suite[J]. Nature protocols, 2016, 11(5): 905-919.

--------- End ---------