几何合理的分片段感知的3D分子生成 FragGen

FragGen 来源于 2024 年 3 月 25 日预印本的文章，文章题目是 Deep Geometry Handling and Fragment-wise Molecular 3D Graph Generation，作者是 Odin Zhang，侯廷军，浙江大学药学院。FragGen 是一个基于分子片段的 3D 分子生成模型。类似的，基于分子片段的 3D 分子生成模型，我们曾经介绍过 Flag。

一、背景介绍

FragGen 来源于浙江大学药学院的崔孙良教授、潘培辰教授、谢昌谕教授和侯廷军教授为通讯作者的文章：《Deep Geometry Handling and Fragment-wise Molecular 3D Graph Generation》文章链接：https://arxiv.org/abs/2404.00014。该文章于 2024 年 4 月 15 日发表在 Arxiv 上。FragGen 的第一作者是 Odin Zhang。我们曾经介绍过他的 3D 分子生成模型 Delete, ResGen。

基于结构的分子图生成方法因为在药物设计中的巨大潜力而备受关注。大多数早期方法采用逐原子生成策略，虽然设计的分子和结构的形状匹配程度较好，结合较好，但往往忽略其他可合成性等关键属性。片段生成策略可以通过组装化学合理的片段来降低合成难度。逐原子生成和片段生成存在的共同问题——它们在协同设计合理的化学和几何结构方面能力有限，导致构象失真和较差的结合倾向。对此，作者引入了深度几何处理协议，将整个几何结构分解成多组几何变量，将设计重点扩展到模型架构之外。

作者讨论了六种协议的优缺点，在此基础上提出了新的混合策略，最终开发了 FragGen ，这是首个可靠的几何片段式分子生成方法。FragGen 在生成几何的质量和分子的合成可行性方面实现了重大突破，解决了分子生成算法应用中的这两个主要挑战。该方法成功生成了纳摩尔级的 II 型激酶抑制剂设计，FragGen 的有效性得到了进一步验证，将其算法进展与现实中的药物开发联系起来。

二、模型介绍

2.1 模型框架

在药物设计中的分子生成模型由于数据稀缺常常导致模型的表现不佳，引入物理约束等先验知识有望解决这一问题。近年来涌现了许多基于结构的分子生成模型，如下图中的 LiGAN、Pkt2Mol、GraphBP、DiffBP 和 ResGen 等。但这些方法在应用时存在一些关键限制：生成的分子构象往往失真，物理上不合理，并且为了适应蛋白口袋常常生成一些包含多个并环的合成性较差的分子。片段式分子生成通过在生成过程中将分子组装成可合成片段，提供了解决多并环问题的方案。现有的成功实现的片段式分子生成方法是 FLAG，但其再每个片段生成步骤的误差会累积，常常导致分子结构的崩溃。因此，一个基于片段的可靠的分子生成方法很重要。

渲染平滑几何是计算物理现实研究的中心焦点，不仅仅是 3D 分子生成，而是几乎所有几何中心应用领域。例如，在分子构象生成中，研究人员采用了“先距离后几何”的协议，首先生成距离矩阵，然后通过在距离约束下优化随机初始化的构象来推导笛卡尔坐标。然而，在将不足以指定的距离矩阵映射到笛卡尔坐标时，非唯一性往往引入额外的误差，导致几何失真。后续的研究探讨了力场优化或端到端的笛卡尔坐标预测，以增强深度学习模型生成准确几何的能力。除了直接生成合理的分子构象外，深度学习在分子对接方面也取得了显著进展。早期模型如 TANKBind 将“先距离后几何”的理念扩展到蛋白质-配体结合构象预测。然而，将蛋白质节点纳入这些模型引入了一个重大挑战：冗余自由度显著增加，导致几何不尽如人意。研究人员随后探讨了端到端的解决方案，直接预测笛卡尔坐标，EquiBind 开创了这一思路。KarmaDock 进一步推动了这一协议，通过采用回收机制模拟经典几何优化，最终将对接成功率提高了约 50% 。然而，所有这些方法仍然在生成不现实的局部结构方面存在困难，如非共面的芳香环和过长的化学键，需要进行后处理步骤，例如几何优化或对齐校正。DiffDock 代表了一种不同的技术方法，专注于调整整体平移、方向和扭转角度等约束变量，以简化分子构象的变换。DiffDock 的理念运作良好，因为它提高了基于深度学习的生成分子的几何合理性，尽管其生成的配体可能仍会与蛋白质口袋残基发生冲突。

在使用深度学习模型正确处理几何问题时，面临的挑战主要有两个方面：几何变量的固有对称性（如下图 A 所示）和几何构建方式。第一方面，对称性考虑，如 SE(3) 不变性/等变性，已经得到了充分解决。许多工作集中在增强模型的特征提取能力，同时强制遵循必要的等变性或不变性原则。第二方面，高层次几何处理协议，例如：EGNN, SchNet, and Geodesic-GNN 等模型。与以对称性为中心的架构设计的发展相比，高层次几何处理协议并未受到同样的关注。计算科学家（在首次进入如药物设计等新领域时）往往倾向于调整模型架构，以期在现有实践下获得更好的表现（例如，给定的几何协议）。但必须认识到，如果目标是实现重大突破，协议本身也应重新评估。所选协议限制了模型的性能边界，并显著影响结果。因此，作者在本项工作中对现有的几何处理协议进行全面审查和重新思考。

作者回顾并总结了六种可用于 3D 分子生成的协议，突出它们各自的挑战，并讨论它们在其他分子几何中心问题中的应用，如分子构象生成和对接问题。在此基础上，提出了一种混合方法，采用多种协议，有效利用每种协议的独特优势，以实现 3D 分子生成的最佳性能，如下图 C 所示。这一新颖策略导致了首个几何可靠且片段式的分子深度生成 FragGen 的开发，如下图 B 所示。实验结果表明，FragGen 达到了最先进的性能，验证了重新制定几何处理协议必要性的论点。此外，作者还将算法开发与现实中的药物设计活动相结合，成功设计了针对白细胞受体酪氨酸激酶的有效 II 型抑制剂（75.9 nM）。这一概念-算法-应用工作不仅填补了基于结构的药物发现中的空白，还丰富了关于几何处理协议的讨论，补充了对称神经网络设计，并为其他与几何相关领域的模型开发提供了蓝图。

2.2 模型性能

作者一共审视了六种通用的几何处理协议，如下图所示，强调了每种协议在考虑蛋白质口袋环境下的三维分子设计中所遇到的独特挑战。内部坐标（Internal Coordinate）协议首先确定三个原子的顺序，然后预测键长、键角和二面角，然而，这种方法往往会导致分子构型的扭曲（逐个原子的生长）。该协议被 GraphBP 方法采用，研究发现其错误主要来源于初始拓扑顺序的错误确定，而在蛋白质口袋环境中正确确定这一顺序本身就非常困难。与 G-SphereNet 等无结构模型不同，在配体单独存在的情况下，拓扑顺序可以自然地沿着生成轨迹展开，而在蛋白质口袋等复杂环境下，内部坐标协议的应用则表现出较大的局限性。相比之下，笛卡尔坐标（Cartesian Coordinate）方法同时生成/预测多个原子的坐标和类型，学习三维坐标概率，但缺乏局部结构约束。这常常导致每个原子位置的错误累积，进而产生不合理的几何形状，例如非共面的环或键长不均匀的苯环。这种挑战在基于扩散模型的方法中尤为明显，如 DiffBP 和 DiffSBDD ，它们通过一次性生成分子，但由于局部约束不足，往往导致几何失真。相对矢量（Relative Vector）协议通过逐个原子生成，预测原子间的坐标向量差异表现得更加稳健。通过确保预测的三维向量满足 SE(3) 等变性，该方法有效地将自由度限制在键长范围内，从而最大限度地减少预测错误对整体几何结构的影响。采用此协议的方法，如 Pocket2Mol 和 ResGen，已经实现了更合理的构象生成。然而，它们仍然面临挑战，特别是在生成多环融合分子时，虽然这些分子有助于增强与蛋白质口袋的相互作用，但其结构复杂且难以合成。

在构像/坐标优化方面，有以下三种方法：

GeomGNN 方法（用于 KarmaDock）利用等变图神经网络学习原子力场，强行更新坐标。此协议因在训练和推理过程简便而受益，因为它避免了不同坐标描述之间的复杂转换。作者在分子生成问题中实现了该方法，结果表明 FragGen-GNN 确实展示了这一优势。然而，它在实现精确原子定位方面仍存在局限性。

GeomOPT 是一种确定下一个原子或片段坐标的经典方法，通过涉及键角和二面角的力场相互作用，理论上可以来避免局部结构的不合理性。尽管该方法有潜力，但它面临显著的限制，包括较长的优化时间以及结构趋向失稳的倾向（陷入局部最优）。

距离几何（Distance Geometry）是另一种用于构象生成的公认方法，如 ConfGF 和 SDEGen 模型中所用，它通过建模原子间距离绕过了神经网络设计中的等变性要求。这降低了模型构建的复杂性，但由于自由度过多，无法从距离矩阵唯一确定三维坐标。因此，即使距离矩阵预测得非常准确，原始笛卡尔坐标的精确重构依然困难，往往导致构象扭曲，如 FLAG 方法中所示。尽管在模型架构设计方面的持续进展旨在提高性能，但它们并未直接解决上述几何协议固有的挑战。认识到这一问题的算法发展不足——这也是影响生成构象整体质量的重要因素之一。

作者提出一种结合策略（Combined Strategy），该策略整合了对每个现有协议优缺点的系统性调查所得的见解。具体来说，结合策略的工作流程如下：作者首先使用相对向量协议进行亚口袋检测，确定后续片段组装的适当位置。在预测下一个片段类型后，其几何结构被分解为局部构象、围绕一个点（连接的原子）的旋转以及围绕一个轴（连接的键）的旋转。传统方法和深度学习方法在处理局部片段几何时通常表现良好。对于围绕一个点的旋转，应用杂化轨道理论约束，例如在标准的 SP3 杂化中（如甲烷中的 109.5° ）一致的键角，以通过基于严格理论见解的化学初始化来指导分子的组装。最后，对于围绕轴的旋转，使用 von Mises 损失直接预测二面角。这种对复杂片段生成几何的解耦方法提供了一种有效的解决方案，后续实验也提供了有力的验证。

2.2.1 FragGen 在 CrossDock 基准测试中的表现

利用新颖的几何处理协议，开发了FragGen，一种基于结构的逐片段分子生成方法。其有效性在广泛认可的 CrossDock 数据集上得到了严格测试，该数据集是先前原子级分子生成研究中的基准。评估包括使用 AutoDock Vina 计算 Vina 评分，以评估配体与靶蛋白的结合亲和力。此外，还包括其他重要的指标，例如药物相似性定量估算（QED）、合成可行性（SA）、Lipinski五规则和辛醇-水分配系数（LogP），以表征生成分子的属性。值得注意的是，SA 成为对比原子级和片段级方法的关键指标，后者通常由于使用现有的商业片段进行组装而产生更高的 SA。基线模型包括四种原子级分子生成方法（GraphBP、DiffBP、Pkt2Mol和ResGen）以及一种片段级模型 FLAG，后者是唯一的开源片段级模型。

每个模型的性能指标如下表所示。FragGen 在 Vina 评分上表现优异，排名如下：FragGen > ResGen > Pkt2Mol > GraphBP > DiffBP > FLAG。FragGen 的 Vina 评分比测试集平均值高出 2.5 kcal/mol，基于热力学原理，转化为超过 100 倍的结合亲和力提升。这个显著的结合效能提升几乎足以将配体的 IC50 从微摩尔级提升到纳摩尔级。此外，FragGen 在生成高质量配体时表现出卓越的化学和几何结构。像 GraphBP 和 DiffBP 这样的原子级方法通常会生成几何结构扭曲的分子，一些由 GraphBP 生成的分子甚至偏离了目标口袋。这些几何结构的缺陷源于内坐标和笛卡尔坐标协议的局限性，其中前者缺乏用于引导生成过程的局部结构约束，后者需要预先定义的拓扑原子顺序。相比之下，使用相对向量协议的 ResGen 和 Pkt2Mol 能够生成更加精确且在视觉上合理的分子几何结构。FLAG 和 FragGen 这两种片段级方法的输出则处于 Vina评分的两端（FLAG: ~-8.9 vs. FragGen: ~-9.9），这是它们几何处理能力的有力证明。FLAG 基于距离几何，常常由于难以将大量成对距离映射到笛卡尔坐标而产生结构不良的分子。相反，FragGen 采用了复杂的几何处理方法，将几何变量分解为四部分，并通过化学知识和端到端学习的结合进行有效管理。

在分子性质方面，FragGen 在 Top 5 结果中获得了 QED 和 SA 的最高评分，强调了其生成分子的化学可行性。这些出色的结果归因于两个关键因素：片段级协议的固有特性和稳健的几何处理方法的优势。片段级协议本质上保证了更好的可合成性，因为它通常将分子分解为一组现有的片段，这也解释了 FLAG 相对较高的 SA 得分。相比之下，像 Pkt2Mol 和 ResGen 这样的原子级方法常常生成完全填充蛋白口袋的分子，导致较低的 QED 和 SA 得分。这种趋势导致药物化学家对将先前的分子生成方法整合到他们的工作流程中持犹豫态度。总之，FragGen 在 Vina 评分、QED 和 SA 方面的进步表明，几何精确性在提高化学合理性方面起着至关重要的作用，因为当前分子状态的几何结构会影响后续片段的结构。对于实际应用来说，FragGen 还确立了自己作为药物发现中有价值工具的地位，尤其是在生成易于合成的样本方面。

2.2.2 FragGen 在著名药物靶标上的表现

为了展示 FragGen 在实际应用场景中的适用性，作者对其在多个著名药物靶标上的表现进行了评估。这些靶标具有良好表征的活性位点和大量实验发现的抑制剂，提供了一个合适的测试平台。与 CrossDock 基准测试不同，本次实验包含了两个额外的分子集合：Active（作为阳性对照的实验验证分子）和Random（从 GEOM-Drug 集合中随机选择的化学片段，作为阴性对照）。下图 A 展示了 FragGen 生成分子（橙色）与基于片段的对照 FLAG（绿色）的结合效力分布。FragGen 的分布更接近于 Active 分子，而 FLAG 则更接近于 Random 集合。这一结果再次强调了片段式分子生成中合理几何协议的优势，准确的几何结构使得与结合蛋白的能量匹配更佳。如下图 B 所示，虽然 ResGen 生成的顶级分子展现出较强的结合效力，但它们在可合成性和药物性方面有所妥协。相反，FragGen 生成的分子不仅在结合效力上与顶级 Active 分子相当（差异仅为约 0.4 kcal/mol），还保持了最高的化学可得性，使其更受化学家的青睐。SA 对比也进一步支持了这一点。

2.2.3 使用 FragGen 设计 LTK 的 II 型抑制剂并进行湿实验验证

激酶是细胞信号传导中的关键酶类，在细胞生长、分化和代谢等多种生理过程中发挥着重要作用。因此，许多激酶抑制剂已经被开发并批准用于治疗癌症、心血管疾病和炎症等疾病。传统的激酶抑制剂（I 型抑制剂）靶向激酶的 ATP 结合位点，虽然具有治疗优势，但在选择性和耐药性问题上存在局限。相比之下，II 型抑制剂（如索拉非尼）靶向另一个别构位点——DFG-out 口袋，可能提供更具选择性且毒性较小的治疗方案。尽管II型抑制剂具有优势，但现有的计算工具（如定量构效关系 QSAR 和对接筛选）在设计超出已知化学空间的有效分子方面存在不足，限制了发现新型治疗药物的可能性。因此，当前的分子生成方法非常适合填补这一空白。

作者选择 LTK 作为验证系统，这是一种有前景的非小细胞肺癌激酶靶点。与以往的回顾性研究不同，通过湿实验进行验证，而不是使用存在争议的对接评分标准。此外，LTK 是一个设计抑制剂较少的新靶点。受到 PDGFRβ 靶点历史药物开发的启发（其基于 I 型抑制剂框架设计 II 型抑制剂），作者使用 FragGen 开发了一个AI驱动的基于结构的工作流程。作者首先基于间变性淋巴瘤激酶（ALK）的蛋白模型构建了 LTK 的 DFG-out 同源模型。然后，将之前报道的 ALK I 型抑制剂对接到 LTK 模型中，旨在通过保留头部铰链结合基团将分子锚定在口袋 I 区域。基于锚定的结构，使用 FragGen 探索针对 II 型口袋的化学空间。在 10 分钟内，FragGen 提出了 97 个化学候选物。随后，应用了四个过滤标准来缩小候选物范围：1）氢供体数量<5；2）氢受体数量<10；3）LogP 在 2 到 5 之间；4）可旋转键数量<10。最终，10 个分子满足这些条件。其中，基于有机化学家的建议，选择了三个用于进一步研究，结构如下图 A 所示。生物测试显示，这些分子对 LTK 具有高亲和力，其中 Darma-1 的效力尤为突出，达到 75.4 nM。其他两个候选物的亲和力分别为 52.4 μM 和 2.56 μM，进一步证明了 FragGen 在蛋白质口袋内设计配体的能力。成功设计出高效的 II 型抑制剂可能归功于 FragGen 对几何结构的精细处理。为说明这一点，作者分析了直接生成的 Darma-1 化合物的结合模式（如下图 B-D 所示）。可以清楚地看到，生成的化合物与 II 型口袋形成了全面的物理相互作用，如与 ASP-155、LYS-35 和 GLU-52 残基形成了三个氢键。如果生成的几何结构不如 FragGen 所提出的合理，分子生成模型即使在对接评分/ADMET评分上表现优异，也会失去实际应用价值：不适当的构象会破坏蛋白质与配体之间的相互作用，降低生成样本的可信度。

2.2.4 生成分子的几何合理性

在 3D 分子生成领域，许多模型依赖几何优化来修正生成分子中的失真，实际上掩盖了深度学习方法在共同设计几何准确分子时的局限性。认识到先前的实验只能间接且定性地解决这些几何挑战，作者引入了两个新指标，以获得更详细和定量的评估：松弛能量（Relax E）和优化后的均方根偏差（OptRMSD）。具体来说，生成的分子经过力场优化后，计算生成分子与优化后分子之间释放的能量和 RMSD ，如下图所示。

下表展示了不同方法的 OptRMSD 和 Relax E 结果。在 OptRMSD 方面，某些模型表现优异。然而，必须承认 OptRMSD 本身对多环结构有偏好。原因在于较大的芳香体系具有更为刚性的框架，不易发生构象变化。因此，诸如 ResGen 和 Pkt2Mol 等倾向于生成多环分子的模型，其较低的 OptRMSD 分数符合预期。相比之下，FragGen 凭借其 OptRMSD 低于 1 Å 的成绩脱颖而出，显示了其生成结构合理分子的能力。

当考虑 Relax E 这一对多环结构偏见较小的指标时，情况有所不同。正如下图 C-D 所示，即使多环结构与简单分子表现出类似的 OptRMSD 值，它们在力场优化后往往会释放更多的能量。在这一背景下，FragGen 再次展示了其优异表现，进一步证明了先前对其几何精度的评估。相反，基于片段的 FLAG 方法以及易于生成失真构象的模型（如 DiffBP 和 GraphBP），在这一指标中表现不佳。

2.2.4 几何处理协议在FragGen中的消融研究

在 3D 分子生成任务中，上述提到的六种协议，有四种——内部坐标、笛卡尔坐标、相对向量和距离几何——已经通过 GraphBP、DiffBP、ResGen 和 FLAG 等工作得到了实现。除了这些，作者还将 GeomGNN 和 GeomOPT 协议整合进 FragGen，从而创建了两个 FragGen 的版本，提供了在 3D 分子生成背景下对每种协议的全面分析。

使用 GeomGNN 协议生成的分子具有最高的结合倾向。然而，这种有利的结合趋势以合成性为代价，相比其他协议降低了约 24% 。这种合成性的降低可归因于模型在填充蛋白质口袋腔体时，局部结构合理性妥协，而没有明确考虑分子的整体合成可行性。另一方面，GeomOPT 方法在合成性方面有显著改善，但在该协议下生成的分子表现出较低的结合倾向。这主要是因为在生成过程中，几何构象被困在蛋白质结构的局部最小值中，导致分子与蛋白质之间的相互作用不佳。综合策略将物理约束与相对向量和内部坐标的优势相结合，成为一种强大的方法。它不仅促进了现实的分子生成，还确保了对目标蛋白质的强结合亲和力。在这种策略下生成的分子不仅表现出更高的结合倾向，超越了所有基线方法（包括原子级和片段级），且在所有协议中展示了最高的合成性。这突显了通过这一综合协议生成的分子结构的有效性和合理性。

三、FragGen 评测

3.1 安装环境

复制代码项目：

git clone https://github.com/HaotianZhangAI4Science/FragGen.git

作者把项目环境压缩到 surfgen.tar.gz 中，上传到了 zenodo 网站，链接是： https://zenodo.org/records/7758282。下载 surfgen.tar.gz，并创建环境，命令如下：

mkdir /workspace/anaconda3/envs/FragGen
tar -xzvf surfgen.tar.gz -C /workspace/anaconda3/envs/FragGen
conda activate FragGen

项目提供训练好的模型，上传到了 zenodo 网站，链接是：https://zenodo.org/records/11119040。一共有三个文件，具体内容如下所示，下载保存到 ./ckpt 中。

项目提供了处理好的 crossdock 的测试集数据，上传到了 zenodo 网站，链接是：Processed files for CrossDock test set 。如下图所示，每个测试体系包含配体结构和蛋白结构以及 ply 格式的信息文件。测试集数据下载后，放在 ./data 文件夹中。

3.2 分子生成案例测试

3.2.1 内置案例

项目的内置案例是 4TOS 蛋白，模型输入需要的处理好的数据保存在 ./data/crosssdock_test/4tos_A_rec_4tos_355_lig_tt_min_0 文件夹中，具体文件如下：

.
|-- 4tos_A_rec.pdb
|-- 4tos_A_rec_4tos_355_lig_tt_min_0.sdf
|-- 4tos_A_rec_4tos_355_lig_tt_min_0_pocket_8.0_res_1.5.ply
`-- info.pkl


0 directories, 4 files

其中，4tos_A_rec_4tos_355_lig_tt_min_0.sdf 和 4tos_A_rec.pdb 分别是配体结构和蛋白结构文件，4tos_A_rec_4tos_355_lig_tt_min_0_pocket_8.0_res_1.5.ply 是预处理的数据信息。

配体在蛋白口袋中的展示如下图：

对应的蛋白表面 ply 文件，可视化如下：

口袋中配体的 2D 结构如下：

作者除了实现组合策略的 FragGen 版本，还把 GeomGNN 和 GeomOPT 协议整合进 FragGen，从而创建了两个 FragGen 的版本，分别是 FragGen-GNN 和 FragGen-OPT。接下来，我们使用这三个版本的模型分别针对内置案例进行分子生成。

3.2.1.1 使用结合策略版本的 FragGen 进行分子生成

使用结合策略版本的 FragGen 进行分子生成的命令如下：

python gen_from_pdb.py \
  --config ./configs/sample_dihedral.yml \
  --surf_file ./data/crosssdock_test/4tos_A_rec_4tos_355_lig_tt_min_0/4tos_A_rec_4tos_355_lig_tt_min_0_pocket_8.0_res_1.5.ply \
  --pdb_file ./data/crosssdock_test/4tos_A_rec_4tos_355_lig_tt_min_0/4tos_A_rec.pdb \
  --sdf_file ./data/crosssdock_test/4tos_A_rec_4tos_355_lig_tt_min_0/4tos_A_rec_4tos_355_lig_tt_min_0.sdf \
  --save_dir ./output_dihedral \
  --device cpu

使用项目提供的配置文件 ./configs/sample_dihedral.yml 。--surf_file 指定基于结构生成的表面信息文件。--pdb_file 和 --sdf_file 分别指定输入的蛋白结构和配体结构。--save_dir 指定生成分子的输出路径。--device 指定使用 cpu 进行采样。

配置文件 ./configs/sample_dihedral.yml 的具体内容如下，设置采样时束搜索的步数，阈值等，基于蛋白口袋采样 100 个分子。

data:
  data_name: 'test'
  dataset:  # test dataset
    name: pl
    path: ./data/crossdocked_pocket10
    split: ./data/split_by_name.pt

model:
    pos_pred_type: dihedral
    checkpoint: ./ckpt/dihedral_208.pt # ./ckpt/freq0_val_90.pt

sample:
  seed: 2020
  mask_init: True
  num_samples: 100 # the maximum number of samples to generate
  beam_size: 500 # beam size is the maximum number of data in the queue, however, with less queue_same_smi_tolorance, it is hard to reach the beam size
  max_steps: 50 # the maximum number of global steps to take
  
  threshold: # the threshold for the first step, low threshold means better diversity, but consumes more time
    focal_threshold: 0.5
    pos_threshold: 0.2
    element_threshold: 0.2
  
  initial_num_steps: 1 # treat the (initial_num_steps) as the initial steps
  next_threshold: # the threshold for following step, low threshold means better diversity, but consumes more time
      focal_threshold: 0.25
      pos_threshold: 0.1
      element_threshold: 0.2
  
  queue_same_smi_tolorance: 2 # high tolerance means better diversity, but consumes more time

设置生成 100 个分子，用时约 4.5 个小时，最终生成 176 个分子，保存在 ./output_dihedral/4tos_A_rec_4tos_355_lig_tt_min_0 文件夹中。

接着，使用 qvina 对接生成分子。创建 ./get_qvina/4ots/FragGen 文件夹，把 ./data/crosssdock_test/4tos_A_rec_4tos_355_lig_tt_min_0/4tos_A_rec.pdb 蛋白结构复制到 ./get_qvina/4ots 文件夹中，把蛋白结构从 pdb 格式转换成 pdbqt 格式，命令如下：

prepare_receptor4.py -r 4tos_A_rec.pdb

创建 ./get_qvina/4ots/FragGen/1_qvina_docking.ipynb 处理生成分子并生成对接脚本。

导入依赖库。创建 ./get_qvina/4ots/FragGen/pdbqt_mol 文件夹，以便保存转换成 pdbqt 格式的生成分子。

import os
import pandas as pd
from rdkit import Chem

if not os.path.exists('./pdbqt_mol'):
    os.makedirs('./pdbqt_mol')

根据参考配体获取对接中心坐标，通过 openbabel 把生成分子转换成 pdbqt 格式的结构，并生成对接脚本 ./get_qvina/4ots/FragGen/docking.sh 。

all_sdfs = [f'/workspace/xxxx/projects/FragGen/output_dihedral/4tos_A_rec_4tos_355_lig_tt_min_0/SDF/{i}.sdf' for i in range(71)]

mol = Chem.SDMolSupplier('/workspace/xxxx/projects/FragGen/data/crosssdock_test/4tos_A_rec_4tos_355_lig_tt_min_0/4tos_A_rec_4tos_355_lig_tt_min_0.sdf')[0]
pos = mol.GetConformer(0).GetPositions()
center = (pos.max(0) + pos.min(0)) / 2

with open('./docking.sh', 'a') as f:

    for i,sdf in enumerate(all_sdfs):

        commands = f"""
    conda activate adt
    obabel {sdf} -O /workspace/xxxx/projects/FragGen/get_qvina/4ots/FragGen/pdbqt_mol/{i}.pdbqt
    """
        os.system(commands)
        f.write(f'''qvina2 \
        --receptor /workspace/xxxx/projects/FragGen/get_qvina/4ots/4tos_A_rec.pdbqt \
        --ligand /workspace/xxxx/projects/FragGen/get_qvina/4ots/FragGen/pdbqt_mol/{i}.pdbqt \
        --cpu 60 \
        --center_x {center[0]} \
        --center_y {center[1]} \
        --center_z {center[2]} \
        --size_x 20 --size_y 20 --size_z 20 \
        --exhaustiveness 16 \
        --log /workspace/xxxx/projects/FragGen/get_qvina/4ots/FragGen/pdbqt_mol/{i}.log\n''')

在终端运行对接脚本，具体命令如下：

conda activate adt
cd /workspace/xxxx/projects/FragGen/get_qvina/4ots/FragGen
sh docking.sh > docking.log

生成分子对接完成后，根据 qvina 的最好对接打分升序排列，保存到 ./get_qvina/4ots/FragGen/Generated_molecules.sdf 中。

import os
from rdkit import Chem

all_logs = [i for i in os.listdir('/workspace/xxxx/projects/FragGen/get_qvina/4ots/FragGen/pdbqt_mol') if i.endswith('.log')]

# 保存生成的分子
mol_sdfs = []

for log in all_logs:
    with open(f'./pdbqt_mol/{log}') as f:
        lines = f.readlines()
    for line in lines:
        if '   1      ' in line:
            aff_num = float([i for i in line.split(' ') if i != ''][1])
            num = log.split('.')[0]
            # 保存生成分子的构象
            mol = Chem.SDMolSupplier(f'/workspace/xxxx/projects/FragGen/output_dihedral/4tos_A_rec_4tos_355_lig_tt_min_0/SDF/{num}.sdf')[0]
            # 给分子添加 vina_score 属性
            qvina_score = aff_num
            mol.SetProp('qvina_score', str(qvina_score))                
            mol_sdfs.append(mol)

# 每个 mol 对象都有 'vina_score' 属性，按 vina_score 降序排列
sorted_mol_list = sorted(mol_sdfs, key=lambda mol: float(mol.GetProp('qvina_score')), reverse=False)

# 把生成分子保存成 sdf 格式
sdf_path = 'Generated_molecules.sdf'
w = Chem.SDWriter(sdf_path)
for mol in sorted_mol_list:
    w.write(mol)
w.close()
print('Generated molecules saved as sdf format in {}!'.format(sdf_path))

所有生成的分子如下：

FragGen_4toz

qvina score 排名前 3 的分子的 2D 结构如下，对应的 qvina score 分数分别为 -12.3， -12.3， -11.9：

qvina score 排名前 3 的分子在口袋中的 Pose 如下，蓝色的是参考分子，紫红色的是生成分子，生成分子相对于参考分子整体上偏小，只占据口袋的部分内部空间。生成的分子与口袋之间的碰撞比较少。

3.2.1.2 使用 GeomGNN 版本的 FragGen 进行分子生成

使用 GeomGNN 版本的 FragGen 进行分子生成的命令如下，使用项目提供的配置文件 ./configs/sample_dihedral.yml 。--surf_file 指定基于结构生成的表面信息文件。--pdb_file 和 --sdf_file 分别指定输入的蛋白结构和配体结构。--save_dir 指定生成分子的输出路径。--device 指定使用 cpu 进行采样。

python gen_from_pdb.py \
  --config ./configs/sample_cartesian.yml \
  --surf_file ./data/crosssdock_test/4tos_A_rec_4tos_355_lig_tt_min_0/4tos_A_rec_4tos_355_lig_tt_min_0_pocket_8.0_res_1.5.ply \
  --pdb_file ./data/crosssdock_test/4tos_A_rec_4tos_355_lig_tt_min_0/4tos_A_rec.pdb \
  --sdf_file ./data/crosssdock_test/4tos_A_rec_4tos_355_lig_tt_min_0/4tos_A_rec_4tos_355_lig_tt_min_0.sdf \
  --save_dir ./output_GeomGNN \
  --device cpu

配置文件 ./configs/sample_cartesian.yml 的具体内容如下，设置采样时束搜索的步数，阈值等，基于蛋白口袋采样 100 个分子。（注：在 FragGen 中不同的生成策略通过配置文件实现）

data:
  data_name: 'test'
  dataset:  # test dataset
    name: pl
    path: ./data/crossdocked_pocket10
    split: ./data/split_by_name.pt

model:
    pos_pred_type: cartesian
    checkpoint: ./ckpt/cartesian_val_158.pt

sample:
  seed: 2020
  mask_init: True
  num_samples: 100 # the maximum number of samples to generate
  beam_size: 500 # beam size is the maximum number of data in the queue, however, with less queue_same_smi_tolorance, it is hard to reach the beam size
  max_steps: 50 # the maximum number of global steps to take
  
  threshold: # the threshold for the first step, low threshold means better diversity, but consumes more time
    focal_threshold: 0.5
    pos_threshold: 0.2
    element_threshold: 0.2
  
  initial_num_steps: 1 # treat the (initial_num_steps) as the initial steps
  next_threshold: # the threshold for following step, low threshold means better diversity, but consumes more time
      focal_threshold: 0.25
      pos_threshold: 0.1
      element_threshold: 0.2
  
  queue_same_smi_tolorance: 2 # high tolerance means better diversity, but consumes more time

设置生成 100 个分子，用时约 42 分钟，最终生成 103 个分子，保存在 ./output_GeomGNN/4tos_A_rec_4tos_355_lig_tt_min_0/SDF 文件夹中。

接着，使用 qvina 对接生成分子。创建 ./get_qvina/4ots/FragGen_GNN 文件夹，创建 ./get_qvina/4ots/FragGen_GNN/1_qvina_docking.ipynb 处理生成分子并生成对接脚本。

导入依赖库。创建 ./get_qvina/4ots/FragGen_GNN/pdbqt_mol 文件夹，以便保存转换成 pdbqt 格式的生成分子。

import os
import pandas as pd
from rdkit import Chem

if not os.path.exists('./pdbqt_mol'):
    os.makedirs('./pdbqt_mol')

根据参考配体获取对接中心坐标，通过 openbabel 把生成分子转换成 pdbqt 格式的结构，并生成对接脚本 ./get_qvina/4ots/FragGen_GNN/docking.sh 。

all_sdfs = [f'/workspace/xxxx/projects/FragGen/output_GeomGNN/4tos_A_rec_4tos_355_lig_tt_min_0/SDF/{i}.sdf' for i in range(103)]

mol = Chem.SDMolSupplier('/workspace/xxxx/projects/FragGen/data/crosssdock_test/4tos_A_rec_4tos_355_lig_tt_min_0/4tos_A_rec_4tos_355_lig_tt_min_0.sdf')[0]
pos = mol.GetConformer(0).GetPositions()
center = (pos.max(0) + pos.min(0)) / 2

with open('./docking.sh', 'a') as f:

    for i,sdf in enumerate(all_sdfs):

        commands = f"""
    conda activate adt
    obabel {sdf} -O /workspace/xxxx/projects/FragGen/get_qvina/4ots/FragGen_GNN/pdbqt_mol/{i}.pdbqt
    """
        os.system(commands)
        f.write(f'''qvina2 \
        --receptor /workspace/xxxx/projects/FragGen/get_qvina/4ots/4tos_A_rec.pdbqt \
        --ligand /workspace/xxxx/projects/FragGen/get_qvina/4ots/FragGen_GNN/pdbqt_mol/{i}.pdbqt \
        --cpu 60 \
        --center_x {center[0]} \
        --center_y {center[1]} \
        --center_z {center[2]} \
        --size_x 20 --size_y 20 --size_z 20 \
        --exhaustiveness 16 \
        --log /workspace/xxxx/projects/FragGen/get_qvina/4ots/FragGen_GNN/pdbqt_mol/{i}.log\n''')

在终端运行对接脚本，具体命令如下：

conda activate adt
cd /workspace/xxxx/projects/FragGen/get_qvina/4ots/FragGen_GNN
sh docking.sh > docking.log

生成分子对接完成后，根据 qvina 的最好对接打分升序排列，保存到 ./get_qvina/4ots/FragGen_GNN/Generated_molecules.sdf 中。

import os
from rdkit import Chem

all_logs = [i for i in os.listdir('/workspace/xxxx/projects/FragGen/get_qvina/4ots/FragGen_GNN/pdbqt_mol') if i.endswith('.log')]

# 保存生成的分子
mol_sdfs = []

for log in all_logs:
    with open(f'./pdbqt_mol/{log}') as f:
        lines = f.readlines()
    for line in lines:
        if '   1      ' in line:
            aff_num = float([i for i in line.split(' ') if i != ''][1])
            num = log.split('.')[0]
            # 保存生成分子的构象
            mol = Chem.SDMolSupplier(f'/workspace/xxxx/projects/FragGen/output_GeomGNN/4tos_A_rec_4tos_355_lig_tt_min_0/SDF/{num}.sdf')[0]
            # 给分子添加 vina_score 属性
            qvina_score = aff_num
            mol.SetProp('qvina_score', str(qvina_score))                
            mol_sdfs.append(mol)

# 每个 mol 对象都有 'vina_score' 属性，按 vina_score 降序排列
sorted_mol_list = sorted(mol_sdfs, key=lambda mol: float(mol.GetProp('qvina_score')), reverse=False)

# 把生成分子保存成 sdf 格式
sdf_path = 'Generated_molecules.sdf'
w = Chem.SDWriter(sdf_path)
for mol in sorted_mol_list:
    w.write(mol)
w.close()
print('Generated molecules saved as sdf format in {}!'.format(sdf_path))

所有生成的分子如下。生成分子相对于参考分子整体上偏小，只占据口袋的部分空间，只在右侧口袋进行了占据。从分子的结构上看，构象确实不如组合策略的那么合理。另外，分子多样性也比较差，可能与参数有关系。与口袋的碰撞分子比较多。

FragGen_GNN_4tos

qvina score 排名前 3 的分子的 2D 结构如下，对应的 qvina score 分数分别为 -12.2， -12.1， -12.0，打分略低于前述的组合策略：

qvina score 排名前 3 的分子在口袋中的 Pose 如下，蓝色的是参考分子，黄色的是生成分子，生成的分子虽然 docking score 很高，但是仅占据了一半的口袋。

3.2.1.3 使用 GeomOPT 版本的 FragGen 进行分子生成

使用 GeomOPT 版本的 FragGen 进行分子生成的命令如下，使用项目提供的配置文件 ./configs/sample_dihedral.yml 。--surf_file 指定基于结构生成的表面信息文件。--pdb_file 和 --sdf_file 分别指定输入的蛋白结构和配体结构。--save_dir 指定生成分子的输出路径。--device 指定使用 cpu 进行采样。

python gen_from_pdb.py \
  --config ./configs/sample_geomopt.yml \
  --surf_file ./data/crosssdock_test/4tos_A_rec_4tos_355_lig_tt_min_0/4tos_A_rec_4tos_355_lig_tt_min_0_pocket_8.0_res_1.5.ply \
  --pdb_file ./data/crosssdock_test/4tos_A_rec_4tos_355_lig_tt_min_0/4tos_A_rec.pdb \
  --sdf_file ./data/crosssdock_test/4tos_A_rec_4tos_355_lig_tt_min_0/4tos_A_rec_4tos_355_lig_tt_min_0.sdf \
  --save_dir ./output_GeomOPT \
  --device cpu

配置文件 ./configs/sample_geomopt.yml 的具体内容如下，设置采样时束搜索的步数，阈值等，基于蛋白口袋采样 100 个分子。

data:
  data_name: 'test'
  dataset:  # test dataset
    name: pl
    path: ./data/crossdocked_pocket10
    split: ./data/split_by_name.pt

model:
    pos_pred_type: geomopt
    checkpoint: ./ckpt/cartesian_val_158.pt

sample:
  seed: 2020
  mask_init: True
  num_samples: 100 # the maximum number of samples to generate
  beam_size: 500 # beam size is the maximum number of data in the queue, however, with less queue_same_smi_tolorance, it is hard to reach the beam size
  max_steps: 50 # the maximum number of global steps to take
  
  threshold: # the threshold for the first step, low threshold means better diversity, but consumes more time
    focal_threshold: 0.5
    pos_threshold: 0.2
    element_threshold: 0.2
  
  initial_num_steps: 1 # treat the (initial_num_steps) as the initial steps
  next_threshold: # the threshold for following step, low threshold means better diversity, but consumes more time
      focal_threshold: 0.25
      pos_threshold: 0.1
      element_threshold: 0.2
  
  queue_same_smi_tolorance: 2 # high tolerance means better diversity, but consumes more time

设置生成 100 个分子，用时约 45 分钟，最终生成 134 个分子，保存在 ./output_GeomOPT/4tos_A_rec_4tos_355_lig_tt_min_0/SDF 文件夹中。

所有生成的分子如下所示：

FragGen_Opt_4tos

接下来，使用类似的 qvina 对接方法，得到 qvina score 排名前 3 的分子的 2D 结构如下，对应的 qvina score 分数分别为 -10.7， -10.5， -10.4（这一打分要低于组合策略以及 GNN 策略）：

qvina score 排名前 3 的分子在口袋中的 Pose 如下，蓝色的是参考分子，橙色的是生成分子，生成分子的质量较差，和蛋白残基有空间冲突的问题，还会生成一些大环结构。

3.2.2 自定义案例

我们使用 3WZE.pdb 蛋白作为自定义的测试案例，模型的输入是参考配体的结构文件和蛋白口袋文件以及基于配体和蛋白结构生成的表面电势文件。蛋白口袋指的是参考配体 8 Å 范围内的残基，使用 pymol 另存为 3wze_pocket.pdb ，配体结构保存为 3wze_ligand.sdf 。创建 ./3WZE 文件夹，把口袋结构 3wze_pocket.pdb 和配体结构 3wze_ligand.sdf 保存到 ./3WZE 文件夹。

配体在口袋中的构象如下：

口袋中分子的 2D 结构，如下：

基于配体和蛋白结构生成的表面电势文件，我们参考 Delete 项目，创建 Surface 环境并使用其中的 surface_maker_test.py 脚本生成表面信息文件。下面介绍 ply 文件的生成方式:

激活 Surface 环境：

conda activate Surface

然后，在项目的主目录 ./ 上运行一下命令：

python surface_maker_test.py \
  --check_software True  \
  --pdb_file ./3WZE/3wze_pocket.pdb \
  --lig_file ./3WZE/3wze_ligand.sdf \
  --outdir ./3WZE/

--pdb_file 指定口袋路径；--lig_file 指定小分子路径；--outdir 制定 ply 文件输出路径。运行结束以后，在 ./3WZE/ 路径下生成 3wze_pocket_pocket_8.0.ply 文件。

至此，模型输入所需的输入数据已经处理完成。接下来使用三个版本的 FragGen 分别在测试案例上生成分子。

3.2.2.1 使用结合策略版本的 FragGen 进行分子生成

所有生成的分子如下：

qvina score 排名前 3 的分子的 2D 结构如下，对应的 qvina score 分数分别为 -11.0， -11.0， -10.2：

qvina score 排名前 3 的分子在口袋中的 Pose 如下，蓝色的是参考分子，紫红色的是生成分子，生成分子和蛋白质残基的空间冲突问题比较严重，只有个别小分子完全位于蛋白口袋中。

3.2.2.2 使用 GeomGNN 版本的 FragGen 进行分子生成

所有生成的分子如下：

FragGen_GNN_3WZE

qvina score 排名前 3 的分子的 2D 结构如下，对应的 qvina score 分数分别为 -9.8， -9.3， -7.8：

qvina score 排名前 3 的分子在口袋中的 Pose 如下，蓝色的是参考分子，紫红色的是生成分子，生成分子相对于参考分子整体上偏小，只占据口袋的部分空间，小分子的构象也非常的差。

3.2.2.3 使用 GeomOPT 版本的 FragGen 进行分子生成

所有生成的分子如下：

FragGen_Opt_3WZE

qvina score 排名前 3 的分子的 2D 结构如下，对应的 qvina score 分数分别为 -10.4， -10.2， -9.9：

qvina score 排名前 3 的分子在口袋中的 Pose 如下，蓝色的是参考分子，紫红色色的是生成分子，生成分子的质量极差，和蛋白残基有空间冲突的问题，还会生成一些并环结构。

通过在三个版本的 FragGen 模型上测试的案例结果来看，使用结合策略的 FragGen 生成的分子相对最好，生成的分子结构中不包含大环，并环等复杂的子结构，但三个版本的生成分子都存在和蛋白质残基空间冲突的问题，这是后续需要改进的地方。

四、总结

在本研究中，作者旨在解决许多三维分子生成模型经常遇到的化学结构和几何结构不合理的问题。通过对现有的六种几何处理协议的分析，提出并开发了 FragGen，这是一种专门为基于结构的片段式分子生成设计的混合策略。

在基准数据集和药物相关靶点上的实验结果表明，FragGen 生成的分子具有最高的结合效力（通过对接评分估算）和合成可行性，满足了实际药物发现工作的需求。几何分析和消融研究表明，FragGen 有效协调了分子几何与蛋白质口袋结构之间复杂的相互作用，凸显了提出的混合策略在结合各种几何处理技术以实现 FragGen 卓越成功中的关键作用。最后，作者成功利用 FragGen 设计了有效的 LTK II 型抑制剂，展示了其实际应用价值，完成了

从概念、算法到应用的完整链条。总而言之，通过整合不同几何处理协议的见解并将其调整为片段式分子生成的特定需求，FragGen 已成为一种稳健的基于结构的药物设计工具。它生成的分子不仅具有高结合亲和力，同时也具备实际的可合成性和几何准确性，这不仅为化学家的工具箱增添了一项实用工具，还作为一种有原则的方法，能够为深度学习中更广泛的几何相关问题推荐合适的几何处理方法。

总体来说，这个是一个很好的概念验证工作，但是距离可以工业上的使用，还有很远的距离。