Nat. Comput. Sci.｜KarmaDock:针对超大规模虚拟筛选的基于深度学习的分子对接方法

本文介绍一篇来自浙江大学侯廷军教授、谢昌谕教授、潘培辰研究员和之江实验室陈广勇研究员团队联合发表的关于分子对接方法的论文。该文章提出了一种基于深度学习的分子对接模型，KarmaDock，可以快速且准确的预测蛋白配体结合构象及其结合强度。该方法通过混合密度函数学习蛋白配体间最优距离分布用于结合强度打分，并将其作为归纳偏置，利用融合自注意力机制的EGNN模型来迭代更新分子坐标从而预测蛋白配体间结合构象。KarmaDock跳过了传统分子对接软件的构象搜索阶段，极大提升了分子对接的速度（0.017 s/complex），适用于超大规模的虚拟筛选。

研究背景

基于结构的的虚拟筛选 (VS) 是药物研发过程中的关键步骤。其中，分子对接技术在预测蛋白质和小分子药物如何结合时起到了关键作用。传统的对接方法，如AutoDock、Glide等，主要通过搜索算法探索配体的可能构象，并用打分函数来评估这些构象。这些方法为了提升虚拟筛选时的效率，对搜索算法和打分函数进行了化简（如采用了启发式搜索算法、打分函数忽略了溶剂效应、引入了经验能量项等），用预测精度换来计算速度。然而，随着分子数据库的日益增大，传统对接方法在超大规模的虚拟筛选面前也显得捉襟见肘。为了解决这个问题，人们通过改良搜索算法（Quick Vina）、通过利用GPU的并行计算能力（AutoDock GPU）或者深度学习算法来提升对接速度。

基于深度学习的分子对接算法主要通过三类方法实现蛋白配体结合构象的预测：（1）预测蛋白配体距离矩阵并重建（如TankBind和EDM-Dock）；（2）类似传统对接算法，利用去噪扩散概率模型预测配体在蛋白质附近的平移、旋转和扭转（DiffDock和FusionDock）；（3）利用EGNN模型预测配体的原子的移动方向和距离并在消息传递过程中更新原子坐标（E3Bind, LigPose和 KarmaDock）。目前，深度学习在保证对接速度和精度的平衡方面仍存在挑战。尽管给定口袋的分子对接在实际应用中更为普遍，但当前主流的深度学习策略主要针对盲对接，因此其应用范围受到了一定限制。大多数模型将结合构象和结合强度的预测视为两个独立的任务，使得在预测蛋白配体结合时无法同时获取结合亲和力，这对于大规模虚拟筛选不够友好。

KarmaDock 架构

KarmaDock的架构包括两个编码器（Graph Transformer (GT) 和 Geometric Vector Perceptrons (GVP)）、一个混合密度网络（MDN）模块进行打分和一个EGNN模块进行对接（见图1 B和方法）。

如图1B所示, 分子被表征为2维的分子图：以原子为节点，共价键为边。蛋白被表征为三维的残基图：以残基为节点并根据距离构建KNN图，K=30。 在进行结合构象预测和打分之前，我们采用了GT编码器和GVP编码器分别编码分子和蛋白内部的相互作用并更新节点特征。

其中H 和 E 分别表示节点和边的特征；S 表示蛋白序列; c, p, s, v 分别表示分子、蛋白、标量和向量。

更新后的节点特征以及节点坐标被用于构建复合物相互做用图, 用于捕捉残基-原子和原子-原子级别的分子间相互作用。在结合构象预测之前，分子的构象由RDKit生成，并在口袋内进行随机旋转和平移。在初始化的时候，我们采用 GraphNorm来初始化相互作用图的节点特征（式3），防止蛋白配体节点的差值过大。然后，节点特征、边特征和分子节点坐标通过堆叠了8层EGNN层的EGNN模块来进行更新并通过recycle3次来增加坐标更新次数（式4）。在每次recycle之前，我们通过门控模块来合并原始节点/边特征和第八层的EGNN输出的节点/边特征（式5）。在EGNN模块结束后，可以选择两种后处理方式（RDKit构象叠合或者力场优化）来保证分子构象的合理性。

F 表示节点/边的特征; p, c, r, l 表示蛋白、分子、 recycle的索引 以及 EGNN 层的索引。

给定节点特征和分子坐标，我们可以根据MDN模块计算结合打分（式6&7）。

其中 H 和 X 分别表示节点特征和节点坐标；p, c, s 分别表示蛋白、分子和标量；

分别表示均值、标准差、混合系数和蛋白配体间距离；score 表示预测的打分值。

KarmaDock的主要创新点如下：(1) 用残基节点来表征蛋白质，以编码几何特征并减少计算成本；(2) MDN模块学习的每个蛋白质和配体节点之间的最小距离的概率分布可以向共享编码器引入距离归纳偏置，从而有助于引导结合构象的学习；(3) 分别设计了两个编码器用于接收距离归纳偏差和学习蛋白质和配体的分子内相互作用；(4) 构建了融合自注意力机制的EGNN的完全连接的交互图以实现快速对接；(5) 采用两种后处理方法，确保生成的构象在键长和键角方面都是合理的。

图1. KarmaDock方法的模型框架

结果与讨论

结合构象预测的精度和速度

作者选择了五种传统对接工具（即AutoDock GPU、Glide SP、Glide XP、Gold和LeDock）进行比较。此外，还在各种测试集上对6个基于深度学习的分子对接模型（即KarmaDock、LigPose、TankBind、E3Bind、EquiBind和DiffDock）的性能进行了比较。值得注意的是，TankBind、E3Bind、EquiBind和DiffDock主要是为盲对接设计的，并使用EquiBind提出的数据集划分方法进行训练和评估。因此，它们的性能不能直接与KarmaDock进行比较。为了确保在EquiBind-split数据集上KarmaDock和这些模型之间的公平比较，作者选择了TankBind作为代表模型，并通过给定真实的结合位点而不是P2Rank预测的结合位点来模拟其在口袋给定情况下的性能。

如表1所示，从精度上来看，KarmaDock的性能在三种数据集划分方法下均超过了传统的对接软件，至少提高了14.9%/22.3%的成功率，与深度学习模型相比，在PDBBind core set和EquiBind划分测试集上精度达到了目前最优（SOTA），在PDBBind refined set 上仅次于LigPose。从速度的角度上看，KarmaDock在PDBBind refined/core集上较传统对接软件实现了至少163.06倍的加速。与其他深度学习模型相比，KarmaDock在速度上仍然有优势，速度至少增加了2.38倍。

此外，为了模拟真实世界虚拟筛选仅有未结合态（APO）蛋白构象可获得的情况，作者还对KarmaDock和传统对接工具在APObind的core set上的性能进行了比较。计算结果表明，KarmaDock在准确性和速度上都超过了所有传统的对接工具。

表1. 传统对接软件和深度学习对接模型的对接精度和速度

结合强度预测的精度

作为一种分子对接工具，KarmaDock不仅要生成结合姿态，还要基于生成的结合姿态给出蛋白质和配体之间的结合强度打分。KarmaDock利用了在RTMScore中已被证明具有令人满意的对接和筛选能力的MDN模块作为打分模块。作者在CASF 2016数据集上测试了KarmaDock作为一个打分函数（SF）的对接能力和筛选能力，这些功能对于虚拟筛选至关重要。计算结果表明，KarmaDock在对接能力和筛选能力方面均表现出强劲的性能，在所有测试的分子对接模型中排名第一。此外，除了RTMScore之外，KarmaDock的表现超过了其他重打分的方法。KarmaDock的性能表明了它强大的构象选择和活性分子富集能力，适用于大规模虚拟筛选。

虚拟筛选能力

在虚拟筛选时，分子对接软件需要准确预测蛋白配体结合构象和结合强度。而深度学习打分函数或没有打分模块的深度学习对接模型无法实现上述目标。KarmaDock作为一种用于分子对接的深度学习模型，在结合构象预测和结合强度预测方面表现出色。本研究评估了KarmaDock与五种传统对接工具（即Glide SP、Gold、LeDock、Surflex和AutoDock Vina）、带有打分模块的深度学习对接模型（TankBind）以及基于序列的模型（TansformerCPI 1.0）在筛选能力评估上的性能。将TransformerCPI纳入评估旨在评估KarmaDock是否可以达到与基于序列的模型相当的速度，且具有高于它们的筛选能力。作者在DEKOIS 2.0上评估了这些模型的性能。从模型角度来看，我们可以从图2中观察到筛选能力排名如下：KarmaDock FF > KarmaDock Raw > KarmaDock Aligned > Glide SP > Surflex > LeDock > Gold ChemPLP > AutoDock Vina > TankBind> TansformerCPI。根据这个排名可以得出四个结论：(1)力场优化可提高KarmaDock 的筛选能力，而RDkit-构象对齐会降低Karma Dock 的筛选能力；(2) Karma Dock 胜过所有测试的传统对接模型以及TankBind 和Tansformer CPI；(3) Tank Bind 和Transformer CPI 排在最后一位，表明预测结合强度而不考虑结合构像无法实现良好的筛选能力；(4) Karma Dock 在筛选能力方面胜过Transformer CPI，速度相似，这表明KarmaDock 可以在超大规模VS 中取代基于序列的模型的作用，通过准确的SFs 或MM/PB(GB)SA 快速过滤可能无效的化合物，并保留其他化合物进行进一步打分。

图2. KarmaDock和其他方法在DEKOIS 2.0 上的富集能力和预测精度

针对LTK靶标的虚拟筛选

在2021年，CLIP1-LTK融合被确认为非小细胞肺癌（NSCLC）的致癌驱动因子，并被视为治疗相应癌症的关键靶点。为了证明KarmaDock在虚拟筛选中的实际意义，我们对LTK靶点进行了真实世界的虚拟筛选。我们使用广泛使用的商业化合物库Specs和ChemDiv数据库进行筛选，这些库包含超过177万个分子。KarmaDock仅用8.4小时就能有效地生成结合构象并打分其结合强度。通过将预测得分最高的配体与晶化ALK PL配体构象之间的相互作用方式进行比较分析，在形成良好相互作用方式的25个配体中选择了一些进行实验验证。在经过测试的25种化合物中，第23号化合物（ChemDiv ID：8005-7327）显示出IC50值为765.6 nM，表明其有潜力成为进一步药物设计的候选化合物。蛋白质和化合物（Lorlatinib和8005-7327）构成的复合物的结合方式以及类似晶体配体化合物在图3 A-C中进行了可视化分析。ALK和Lorlatinib（图3 A）的结合构象分析发现，Lorlatinib与GLU-1197、MET-1199和ARG-1120形成氢键，并与GLY-1269和ASP-1270形成卤素键。KarmaDock生成的LTK与Lorlatinib之间的结合也呈现类似相互作用模式（图3 B），即与GLU-591和MET-593形成氢键，并与GLY-663和ASN-648形成卤素键，这表明KarmaDock能够基于相互作用准确地生成PL的结合位姿。如图3 C所示，选择的第23号化合物（ChemDiv ID：8005-7327）也与MET-593形成氢键，类似于前两个复合物，这有助于抑制LTK并暗示KarmaDock认为LTK与化合物之间的氢键是必不可少的。此外，在该化合物上存在芳香性氢键，在提高结合强度方面起到了进一步作用。可以基于该化合物的结构进行进一步设计和优化，以探索其作为药物候选化合物的潜力。

图3. 针对LTK靶标的虚拟筛选和实验验证：（A）靶标ALK和化合物Lorlatinib在PDB中的结合构象和相互作用模式（PDB ID：4CTB）；（B）由KarmaDock预测的LTK和化合物Lorlatinib的结合构象和相互作用模式（LTK的蛋白质结构由Alphafold 2生成）；（C）由KarmaDock预测的靶标LTK和化合物23（ChemDiv ID：8005-7327）的结合位姿和相互作用模式；(D)(E)分别是Lorlatinib和化合物23 (ChemDiv ID: 8005-7327) 的吸收曲线。Ba/F3-CLIP1-LTK细胞抑制活性显示为三个独立实验的均值±标准偏差。

总结

作者提出了一种名为KarmaDock的深度学习模型，能够快速、精确的生成结合构象并预测结合强度。并且通过后处理技术，如力场（FF）优化和RDkit构象对齐，以纠正预测配体构象中检测到的键长和角度不准确问题。

值得注意的是，这些后处理方法确实提高了构象合理性，但却降低了对接精度，在两者之间取得平衡成为一个具有挑战性的任务。一个务实的解决方案是对化合物库进行分层筛选，在筛选过程中首先忽略构象合理性，并随后对得分最高的复合物进行FF优化。进一步进行深入实验验证可以确认这个解决方案，并扩大KarmaDock应用范围。

然而，作者也指出KarmaDock存在固有局限性，它作为一种半柔性对接工具在蛋白质结构变异方面表现出盲目性。这个特点或许可以解释APObind上观察到的性能下降。

参考资料

Zhang, X.; Zhang, O.; Shen, C.; Qu, W.; Chen, S.; Cao, H.; Kang, Y.; Wang, Z.; Wang, E.; Zhang, J.; Deng, Y.; Liu, F.; Wang, T.; Du, H.; Wang, L.; Pan, P.; Chen, G.; Hsieh, C.-Y.; Hou, T., Efficient and accurate large library ligand docking with KarmaDock. Nature Computational Science 2023, 3, 789-804

https://www.nature.com/articles/s43588-023-00511-5