GCDM 是一个新的 3D 分子生成扩散模型，与之前的 EDM 相比，GCDM 优化了其中的图神神经网络部分，使用手性敏感的 SE3 等变神经网络 GCPNET++ 代替了 EDM 中的 EGNN，让节点间消息传递、聚合根据手性不同而进行。本文对 GCDM-SBDD，即口袋条件下的分子生成进行了测试。

一、背景介绍

GCDM 来源于美国密苏里大学电气工程与计算机科学系的研究助理 Alex Morehead 为通讯作者的文章：《Geometry-complete diffusion for 3D molecule generation and optimization》。文章链接：https://www.nature.com/articles/s42004-024-01233-z 。该文章于 2024 年 7 月 3 日发表在 《 communications chemistry 》上。

近来，扩散模型和等变图神经网络（GNNs）经常被用于生成三维分子。然而，这些方法无法学习三维分子的重要几何性质，因为它们采用了与分子无关的非几何图神经网络（GNN）作为其 3D 图去噪网络，这明显阻碍了它们生成有效的大型 3D 分子。在本研究中，作者通过引入用于 3D 分子生成的几何-完备扩散模型（Geometry-Complete Diffusion Model， GCDM）来解决这些问题。GCDM 在条件和非条件设置下，分别在 QM9 数据集和更大的GEOM-Drugs 数据集上表现出了比现有 3D 分子扩散模型显著更好的性能。重要的是，作者展示了 GCDM 的生成去噪过程使模型能够生成大量有效且能量稳定的大型分子，而现有方法在所学习的特征上未能实现这一点。此外，还展示了 GCDM 的扩展不仅能够有效设计针对特定蛋白质口袋的 3D 分子，还能够被重新用于一致地优化现有 3D 分子的几何结构和化学组成，以提高分子的稳定性和特性特异性，展示了分子扩散模型的新多功能性。

二、模型介绍

近年来，扩散概率模型 (Diffusion Probabilistic Models， DPMs) 在图像生成、音频合成等领域取得了巨大成功。同时，几何深度学习在蛋白质结构预测等科学问题上也展现出强大潜力。然而，将这两种方法有效结合以解决 3D 分子生成的挑战仍然是一个问题。

现有的 3D 分子生成无法有效学习 3D 分子的重要几何特性，难以生成有效的大型 3D 分子结构，缺乏对分子手性的考虑，在条件生成和分子优化任务上表现不佳。为了克服这些限制，作者提出了用于 3D 分子生成的几何-完备扩散模型（Geometry-Complete Diffusion Model，GCDM），通过引入几何完备的图神经网络和手性敏感的去噪过程，实现高质量 3D 分子的生成和优化。

GCDM 的核心思想是将几何完备的图神经网络与扩散概率模型相结合，实现对 3D 分子结构的精确建模和生成。具体来说，GCDM 具有以下关键特点：

（1）几何完备性：GCDM 采用 GCPNET++ 作为骨干网络，能够通过手性敏感的消息图同时处理标量(不变)特征和向量(等变)特征，从而捕获 3D 分子的完整几何信息。

（2）手性敏感: 模型通过引入手性敏感的局部坐标系，能够准确建模分子的手性特性。

（3）SE(3) 等变性: GCDM 保持了对 3D 旋转和平移变换的等变性，确保生成的分子结构在空间变换下保持一致。

（4）联合扩散过程: 模型对原子坐标和原子类型进行联合扩散，实现了对分子拓扑和几何结构的统一建模。

（5）条件生成能力: GCDM 支持基于分子性质或蛋白质口袋的条件生成，使其能够应用于更广泛的实际场景。

在本研究中，作者通过一系列实验评估结果得到了如下结论：

（1）通过几何执行消息传递的神经网络使 3D 分子的扩散生成模型能够生成有效且能量稳定的大分子，而非几何消息传递网络无法做到这一点，作者还引入了关键计算指标来支持这些发现。

（2）物理归纳偏差（如不变图注意力和分子手性）在生成有效 3D 分子方面都起着重要作用。

（3）新提出的几何完备扩散模型（GCDM），是首个结合上述见解并实现 3D 分子生成所需的理想等变性（即SE(3)等变性）的扩散模型，在 QM9 数据集的条件 3D 分子生成和 GEOM-Drugs 数据集的大 3D 分子非条件生成方面均创下了最新的最先进（SOTA）成果。在后者中，GCDM 的 PoseBusters 有效率超过两倍；在 QM9 数据集上在非条件场景下生成了更多独特且新颖的小分子；并在蛋白质条件 3D 分子生成中实现了更好的 Vina 能量分数和超过两倍的 PoseBusters 有效率。

（4）作者进一步证明，几何扩散模型如 GCDM 可以持续进行分子的 3D 优化，无需重新训练即可优化分子的稳定性和特定分子特性，而非几何扩散模型则无法做到。

2.1 模型框架

下图为几何完备性扩散模型 (GCDM) 的整体框架。GCDM 引入了改进版的 SE(3)-等变图卷积网络（GCPNET++），解决了现有方法在处理 3D 分子几何特性时的不足。GCPNET++ 能够处理标量（不变）和向量（等变）节点和边特征，利用消息传递机制在图神经网络中，通过节点和边之间的信息交换来学习图的结构和属性，其中，图中的消息传递考虑手性。具体来说，几何量传递（如原子间距离、角度等）使得模型能够学习并保持分子的几何特性。此外，模型还结合了物理归纳偏置，包括不变图注意力和分子手性，以确保生成分子的物理化学属性与真实分子一致。

注：GCPNET++ 中的图不是一个简单的 EDM 图（不是简单基于距离构建的图）

在去噪扩散过程中，GCDM 包括两个主要步骤：正向过程和反向过程。在正向过程中，逐步向输入分子图添加噪声，使得分子图逐步退化为纯噪声。而在反向过程中，通过训练模型学习去噪过程，从噪声状态逐步恢复出原始分子图。具体步骤包括从噪声数据分布采样，并通过迭代去除模型预测的噪声来生成新的分子。这一过程中，使用了零重心技巧和等变神经网络，确保模型生成的分子具有平移和旋转不变性。

训练过程旨在最小化负对数似然（NLL），以优化模型的生成能力。首先，对 QM9 和 GEOM-Drugs 数据集进行预处理，包括原子坐标归一化、分子图的构建等。然后，通过噪声参数化在模型输入上添加固定的 Markov 链方差调度噪声。损失函数用于最小化预测噪声与实际噪声之间的差异，从而优化模型参数。在训练过程中，使用等变神经网络确保生成分子的平移和旋转不变性。此外，还通过物理归纳偏置，确保生成分子的物理化学属性与真实分子一致。

GCDM 的技术实现主要包括以下几个关键组成部分：

（1）GCDM 将 3D 分子使用全连接的图表示，包含节点（原子）信息和节点（原子）的连接情况。

（2）GCDM 定义了一个联合扩散过程，同时也对等变的原子坐标和不变的原子类型进行扩散。

（3）GCDM 是通过预测噪音来参数化逆过程的。

（4）GCDM 采用 GCPNET++ 作为消噪网络，该网络能够同时学习标量(不变)和向量(等变)节点和边特征，通过手性敏感的图消息传递过程实现几何完备的表示学习。

（5）基于先前扩散模型，GCDM 噪音参数化后得到的优化目标为：

以实现预测噪音尽可能地逼近真实噪音。

2.2 评估指标、数据集和基线模型

作者通过计算在相应测试数据集上的平均负对数似然（NLL），来反映预测原子类型和坐标的去噪配对效果，数值越低表明预测效果越好。在分子评估层面，利用原子对之间的距离及其各自的原子类型来预测键类型（单键、双键、三键或无键），测量了生成的原子中具有正确化合价的比例（原子稳定性，AS）和生成的所有原子均稳定的分子的比例（分子稳定性，MS）。为了进一步评估方法在 3D 分子生成方面的表现，作者还通过 RDKit 计算分子有效性（Val）和唯一性（Val and Uniq）以及通过 PoseBusters 测试的 PoseBusters 有效性（PB-Valid），PB-Valid 包括各类化学和结构有效性测试，如所有原子连接、键长和键角有效、无内部空间冲突、平面芳香环和双键、低内部能量、正确化合价和可 kekulize 化等。

实验使用的 QM9 数据集包含 130k 个小分子的分子属性和 3D 原子坐标。在 QM9 中的每个分子最多包含 29 个原子，数据集预处理后为每个分子补充氢原子。对于 3D 分子生成任务，训练GCDM 来无条件生成分子，生成每个分子的原子类型（H、C、N、O 和 F）、整数原子电荷以及这些原子的 3D 坐标。将 QM9 数据集划分为训练集、验证集和测试集，分别包含 100k、18k 和 13k 个分子示例。

实验中使用的 GEOM-Drugs 是一个知名的大型 3D 分子构象数据集，用于下游机器学习任务。该数据集包含 43 万分子，每个分子平均有 44 个原子，经过数据集处理后，每个分子最多可以包含 181 个原子。在这次实验中，作者收集了每个分子对应的 30 个最低能量的构象，并要求每个基线方法生成这些分子的 3D 位置和组成原子的类型。

作者选择了 10 种用于 3D 分子生成的基线模型和 GCDM 进行了比较，每个模型均在相同的QM9 划分数据上训练和测试，以确保公平比较：包括G-Schnet 、等变归一化流（E-NF） 、图扩散模型（GDM） 及其变体（即 GCM-aug）、等变扩散模型（EDM） 、Bridge 和 Bridge+Force 、隐向量扩散模型（LDMs），如 GraphLDM 及其变体 GraphLDM-aug ，以及最先进的 GeoLDM 方法 。作者选择这些基线模型作为隐式键预测方法的代表，这些方法通过生成分子的原子类型和原子间距离推断键类型。为了更深入分析 GCDM 某些关键模型组件的影响，作者还包括了两个 GCDM 的消融模型。这两个消融模型分别为没有手性和几何完备局部框架的 GCDM（即 GCDM w/o Frames）以及没有对每条边消息应用标量消息注意力（SMA）的GCDM（即 GCDM w/o SMA）。

2.3 模型性能

作者在多个数据集和任务上评估了 GCDM 的性能，包括 QM9、GEOM-Drugs、Binding MOAD 和 CrossDocked 数据集。下面按照不同数据集和任务场景下分别介绍评估结果。

2.3.1 QM9 数据集上的无条件 3D 分子生成

下图是 GCDM 和基线模型的 3D 分子生成评估结果。在表的上半部分，可以看到 GCDM 相比所有基准方法在可能性（NLL）、有效性和唯一性方面的分子生成表现均优于所有基线方法，虽然 AS 和 MS 结果略低于 GeoLDM，但标准偏差较低。在表的下半部分，作者使用 5 次采样运行重新评估 GCDM 和 GeoLDM，并报告每个指标的 95% 置信区间。

与 GeoLDM 相比，GCDM 生成了 1.6% 更多的 RDKit 有效且唯一的分子，并且生成了多 5.2% 的新颖分子，同时在 QM9 测试数据集上表现出最佳的 NLL。这一结果表明，虽然 GeoLDM 在新颖性率上接近平衡（即50%），但 GCDM 在平均生成的分子中几乎 60% 的情况下能接近 GeoLDM 的稳定性和 PB 有效率，表明 GCDM 可能在准确探索新颖且有效的小分子空间方面更有用。

GCDM 的 SMA 消融实验表明，为生成稳定的 3D 分子，GCDM 高度依赖于能够执行一种轻量级的完连接图自注意力机制，这表明未来研究需要将这种生成模型扩展到蛋白质等大型生物分子。此外，从 GCDM 中移除几何局部框架嵌入揭示了分子手性和几何完备性的归纳偏差是 GCDM 取得这些 SOTA（最先进）结果的重要贡献因素。

下图展示了 GCDM 生成的 QM9 大小分子的 PoseBusters 有效示例。

2.3.2 QM9 数据集上的条件生成

为了实现 3D 分子条件生成的实际应用，作者将 GCDM 与现有的 E(3) 等变模型 EDM 和GeoLDM 进行比较，还与两个简单的基线进行比较：“Naive (Upper-bound) ”，即通过分子属性分类器 ϕ_{c}根据生成方法生成的 3D 分子和打乱（即随机）的属性标签来预测分子属性；以及“# Atoms”，即通过生成方法生成的 3D 分子的原子数来预测其分子属性。对于每种基线方法，报告了由三个 EGNN 分类器 ϕ_{c}组成的集成模型在分子属性预测方面的平均绝对误差（MAE）。对于 GCDM，通过将其条件模型训练在六种不同的分子属性特征输入（α、gap、homo、lumo、μ 和 C_{v}）上，使用 QM9 验证集作为模型的训练数据集，并使用 QM9 训练集作为相应 EGNN 分类器集成的训练数据集进行训练。因此，可以预期当一种方法生成的属性特定分子越准确时，其性能与“QM9 (Lower-bound)”的差距就会越小。

从下表可以看出，与所有基线模型相比，GCDM 在给定分子属性的条件生成方面取得了最佳的整体结果。特别是表格下半部分，GCDM在所有六种分子属性的 MAE 结果上平均超过了 SOTA GeoLDM 方法，分别为28%、9%、3%、15%、21%和35%，同时几乎和 GeoLDM 的 PB-Valid 率相当。这些结果定性和定量地表明，使用几何完备扩散，GCDM能够显著精确地生成具有特定分子属性（例如α-极化率）的3D分子。

条件生成的样本如下图所示，分子是通过 GCDM 生成的，并使用不断增加的 α 值。随着 α 从 68.7 (a) 到 93.6 (f) 变化，生成分子的结构特征逐渐发生改变。

2.3.3 GEOM-Drugs 数据集上的大分子生成

首先，下表显示了一个值得注意的现象：由于 GEOM-Drugs 分子的尺寸和原子复杂性，基于这些原子间距离估算键类型所积累的误差，导致基线方法在此数据集测量的分子稳定性指标比 QM9 数据集低得多。因此，为了准确评估方法在这种情况下的性能，报告额外的化学和结构有效性指标（例如PB-Valid）至关重要，作者在表下半部分做了比较。GCDM 在 GEOM-Drugs 数据集上超越了 EDM 的 SOTA 负对数似然结果，提升了57%，并将 GeoLDM 的 SOTA 原子和分子稳定性结果分别提升了 4% 和六倍以上。更重要的是，GCDM 能够生成大量 PB 有效的大分子，甚至超越了 GEOM-Drugs 数据集的参考分子稳定性率，提升了 54%，表明几何扩散模型（如 GCDM）不仅能够有效生成有效的大分子，还可以超越 GEOM-Drugs 中稳定分子的本地分布进行推广。

下图展示了 GCDM 在 GEOM-Drugs 上生成的 PB 有效的大分子示例。作为 GCDM 为生成的分子图产生低能量结构概念的一个例子，图中 a 和 f 的自由能分别使用 GFN2-XTB 理论水平下的 CREST 2.1239 计算为 -3 kcal/mol 和 -2 kcal/mol ，这与 GEOM-Drugs 数据集的对应自由能分布均值（-2.5 kcal/mol）相匹配。

为检测一种方法在总体上是否生成具有不太可能的 3D 构象的分子，生成分子的能量比定义为分子 UFF 计算能量和同一分子图的 50 个 RDKit ETKDGv3 生成构象的平均值的比率。生成分子能量比大于 7 被认为具有高度不可能的 3D 构象。下图显示了 GCDM 生成的大型 3D 分子的平均能量比明显更低且边界更紧密，与此任务的最新基线模型 GeoLDM 相比，GCDM 还生成了更具能量稳定性的 3D 分子构象。

2.3.4 QM9 数据集上的分子性质引导的 3D 分子优化

为了评估分子扩散模型是否不仅能够生成新的 3D 分子，还可以在分子属性引导下优化现有的小分子，作者采用了 QM9 数据集进行实验。首先使用无条件 GCDM 模型通过 10个 时间步逆扩散生成 1000 个 3D 分子（使这些分子处于未优化状态），然后将这些分子提供给一个独立的属性条件扩散模型，以优化这些分子的原子类型和 3D 坐标。该条件模型接受这些 3D 分子作为中间状态，并在 100 和 250 个时间步内进行属性引导优化。最后使用一个外部属性分类器对这些优化后的分子进行评分，评估（1）优化后的分子预测的属性值是否得到了相应的改善（第一个指标），以及（2）优化过程中分子的稳定性是否发生了变化（第二个指标）。

作者定量探讨生成模型是否能够降低属性特异性 MAE 并提高现有 3D 分子稳定性的效果，下图中展示了一个实际发现：几何扩散模型（如 GCDM）可以通过最小的修改重新用作 3D 分子优化方法，提高分子的稳定性和属性特异性。这一发现实验证明，分子去噪扩散模型可以在药物发现流程的优化阶段中应用，以比以往更快的速度实验更广泛的潜在药物候选（后优化）。

同时，基线 EDM 方法未能一致地优化现有 3D 分子的稳定性和属性特异性，这表明几何方法如 GCDM 在理论和经验上更适合此类任务。平均而言，GCDM 在 100 步时间尺度上将初始分子的稳定性提高了超过 25% ，将每个分子属性的特异性提高了超过 27% ，而 EDM 在 100 步时间尺度上将分子的稳定性提高了 13% ，将属性特异性提高了 15% 。将优化时间步数从 100 步增加到 250 步不一致地进一步改善了分子的稳定性和属性特异性，这表明优化轨迹可能在 100 步左右达到局部最小值，因此合理化了将 1000 个分子的优化计算时间从 15 分钟（250步）减少到 5 分钟（100步）的操作。

2.3.5 蛋白质条件下的3D分子生成

为了研究几何完整方法是否能够增强分子扩散模型在给定蛋白质口袋中生成 3D 模型的能力（即进行基于结构的药物设计（SBDD）），在本次实验中，作者采用了标准的 Binding MOAD 和 CrossDocked 数据集对几何完整扩散生成模型 GCDM-SBDD 进行训练和评估。GCDM-SBDD 基于 GCPNET++，用于蛋白质口袋感知分子生成。Binding MOAD 数据集包含100000 个高质量的蛋白质-配体复合物用于训练，并有 130 个蛋白质用于测试，使用 30% 的序列同一性阈值定义交叉验证的划分。同样，CrossDocked 数据集包含 40484 个高质量的蛋白质-配体复合物，划分为训练集（40354）和测试集（100）。

在该任务上，选择 DiffSBDD-cond 和 DiffSBDD-joint 作为基线模型，与几何完备、蛋白质感知扩散模型 GCDM-SBDD 进行对比。使用 QuickVina 打分、药物相似性（QED）、合成可及性（SA）、PB-Valid、平均满足 Lipinski 五规则的条数以及分子平均多样性作为评价指标。“joint” 和 “cond” 配置分别代表针对蛋白质目标生成分子时，是否同时修改蛋白质目标内结合口袋的坐标。

下图显示，在两个标准的 SBDD 数据集（即 Binding MOAD，简称 BM 和 CrossDocked，简称 CD）上，GCDM-SBDD 生成的分子相比现有方法具有更少的冲突（PB-Valid）且能量更低（Vina）。此外，在所有其他指标上，GCDM-SBDD 在药物相似性指标（如 QED）方面取得了相当或更好的结果。这些结果表明，GCDM 结合 GCPNET++ 作为其去噪神经网络，不仅在全新 3D 分子生成方面表现良好，还能针对蛋白质目标生成特定的 3D 分子，显著扩展了 GCDM 的实际应用领域。具体而言，GCDM-SBDD 在两个数据集上将 DiffSBDD 的平均 Vina 能量分数提高了 8%，并在更具挑战性的 Binding MOAD 数据集上生成了超过两倍的 PB 有效“候选”分子。

PB 有效比例在考虑和不考虑蛋白质-配体空间冲突时的差异表明，基于深度学习的药物设计方法在生成针对特定蛋白质口袋的分子后，可能会显著受益于交互感知的分子动力学松弛，这可能允许许多生成的“候选”分子通过这种松弛恢复 PB 有效性。

下图展示了 GCDM 能够为未见过的蛋白质目标生成无碰撞的真实且多样化的低 Vina 能量的 3D 分子。通过 GCDM-SBDD 生成的蛋白质口袋特异性 3D 分子。GCDM-SBDD 为Binding MOAD (a， b) 和 Cross Docked (c， d) 测试蛋白生成的分子。这些选择的口袋结合分子的 Vina 能量评分范围从 −8.2 (c) 到 −9.7 (b)。

三、GCDM 模型评测

文章中，作者提供了 GitHub 的链接，但是这个链接是针对纯的分子生成的。为了测评在药物设计的能力，我们选择的是针对口袋生成的 GitHub 代码。

3.1 安装环境

复制代码项目：

git clone https://github.com/BioinfoMachineLearning/GCDM-SBDD.git

通过项目提供的 env.yml 创建 GCDM-SBDD 环境

conda env create -f environment.yaml
conda activate GCDM-SBDD
#  install local project as package
pip3 install -e .

在环境中安装 QuickVina 2 ，以便后续对接。

wget https://github.com/QVina/qvina/raw/master/bin/qvina2.1
chmod +x qvina2.1
# 将 qvina 迁移到 conda 环境内
mv qvina2.1 /anaconda3/envs/GCDM-SBDD/bin

下载模型的 checkpoint 文件夹

# fetch and extract model checkpoints directory
wget https://zenodo.org/record/13375913/files/GCDM_SBDD_Checkpoints.tar.gz
tar -xzf GCDM_SBDD_Checkpoints.tar.gz
rm GCDM_SBDD_Checkpoints.tar.gz

./checkpoint 文件夹中包含已经训练好的模型，具体内容如下：

.
|-- binding_moad_ca_joint_egnn.ckpt
|-- binding_moad_ca_joint_gcpnet.ckpt
|-- bindingmoad_ca_cond_egnn.ckpt
|-- bindingmoad_ca_cond_gcpnet.ckpt
|-- crossdocked_ca_cond_egnn.ckpt
|-- crossdocked_ca_cond_gcpnet.ckpt
|-- crossdocked_ca_joint_egnn.ckpt
`-- crossdocked_ca_joint_gcpnet.ckpt


0 directories， 8 files

一共包括 8 个 checkpoint 文件。

3.2 分子生成案例测试

3.2.1 内置案例

项目提供的基于口袋的从头分子生成案例使用的是 4OZ2 蛋白，原始配体在口袋中的结构如下(原始配体小分子名字为 1Z6 ，A:501)：

配体 1Z6 的 2D 结构如下：

./notebooks 文件夹中提供了两个 jupyter notebook （即 pocket_based_molecule_generation.ipynb 和 molecule_filtering_with_posebusters.ipynb）分别指导（1）基于口袋的分子生成和（2）使用 PoseBusters 组件过滤生成分子。接下来，按照项目提供的教程对内置案例进行分子生成和测试。

3.2.1.1 分子生成

首先，根据 pocket_based_molecule_generation.ipynb 的教程准备数据，生成分子。通过下面命令，从 PDB 数据库下载 4OZ2.pdb 蛋白结构

wget https://files.rcsb.org/download/4OZ2.pdb1.gz && gunzip 4OZ2.pdb1.gz && mv 4OZ2.pdb1 4OZ2.pdb

创建 ./experiments 文件夹，在其中继续创建一个根据蛋白配体复合物命令的子文件夹，即 ./experiments/4OZ2_1Z6_07_25_2024_12_00_00 。案例的下载文件 4OZ2.pdb  放在该文件夹中。创建案例的输出文件夹 experiments/4OZ2_1Z6_07_25_2024_12_00_00/outputs_07_25_2024_12_00_00/，将用于保存生成的分子，创建命令如下：

mkdir -p experiments/4OZ2_1Z6_07_25_2024_12_00_00/outputs_07_25_2024_12_00_00/

通过 MOE 把配体-蛋白复合物结构 4OZ2.pdb 中的参考配体 1Z6 删除，保存为 4OZ2_without_1Z6.pdb。4OZ2.pdb 和 4OZ2_without_1Z6.pdb 两个结构文件放在 ./experiments/4OZ2_1Z6_07_25_2024_12_00_00 文件夹中，供后续评价生成和评价分子使用。

基于 4OZ2 的蛋白口袋生成一批适应配体 1Z6 结合位点的 3D 分子。具体命令如下：

python3 generate_ligands.py \
  checkpoints/bindingmoad_ca_cond_gcpnet.ckpt \
  --pdbfile experiments/4OZ2_1Z6_07_25_2024_12_00_00/4OZ2.pdb \
  --outdir experiments/4OZ2_1Z6_07_25_2024_12_00_00/outputs_07_25_2024_12_00_00/ \
  --ref_ligand A:501 \
  --n_samples 100 \
  --num_nodes_lig 13 \
  --sanitize \
  --relax

指定模型为 checkpoints/bindingmoad_ca_cond_gcpnet.ckpt ；

--pdbfile experiments/4OZ2_1Z6_07_25_2024_12_00_00/4OZ2.pdb 指定下载的蛋白质-配体复合物位置；

--outdir experiments/4OZ2_1Z6_07_25_2024_12_00_00/outputs_07_25_2024_12_00_00/ 指定生成分子的保存路径；

--ref_ligand A:501 指定参考配体 1Z6 的信息，包括所属的 A 链，index 为 501 ；

--n_samples 100 指定采样生成 100 个分子；

--num_nodes_lig 13 指定参考配体的节点数目，即非氢重原子数目；

--sanitize 表示确保生成分子的结构有效性，检查分子是否符合化学规则，包括环的闭合、化学价键等；

--relax 表示使用 UFF 力场优化生成分子 200 步。

运行输出：

alphas2 [9.99990000e-01 9.99988000e-01 9.99982000e-01 ... 2.59676966e-05
 1.39959211e-05 1.00039959e-05]
gamma [-11.51291546 -11.33059532 -10.92513058 ...  10.55863126  11.17673063
  11.51251595]
Entropy of n_nodes: H[N] 7.115581512451172
... ...

生成分子保存为 ./experiments/4OZ2_1Z6_07_25_2024_12_00_00/outputs_07_25_2024_12_00_00/4OZ2_mol.sdf ，设置生成 100 个分子，实际生成 84 个分子。num_nodes_lig 和参考配体 1Z6 的非氢重原子数目一样，为 13 。生成过程显存占用约 10 GB，用时约 20 分钟。

3.2.1.2 通过 PoseBusters 评估生成分子

根据 molecule_filtering_with_posebusters.ipynb 的教程评价并过滤分子。创建 ./experiments/4OZ2_1Z6_07_25_2024_12_00_00/pb_screen_mols.py 脚本，通过 PoseBusters 筛选生成分子。下面简要介绍脚本的代码。

导入筛选分子所需的依赖库。

import argparse
import glob
import os
import shutil
import tempfile


import pandas as pd


from posebusters import PoseBusters
from rdkit import Chem

选定 PoseBusters 组件筛选分子的测试类型。包括化学有效性（分子是否是化学上合法的）、分子内有效性（检查键长键角、空间冲突、内部能量等）、分子间有效性（分子与蛋白、辅因子的距离是否合适，是否有空间冲突）和 蛋白-配体的最大距离等方面。

BUST_TEST_COLUMNS = [
    # chemical validity and consistency #
    "mol_pred_loaded",
    "mol_cond_loaded",
    "sanitization",
    # intramolecular validity #
    "all_atoms_connected",
    "bond_lengths",
    "bond_angles",
    "internal_steric_clash",
    "aromatic_ring_flatness",
    "double_bond_flatness",
    "internal_energy",
    # intermolecular validity #
    "minimum_distance_to_protein",
    "minimum_distance_to_organic_cofactors",
    "minimum_distance_to_inorganic_cofactors",
    "volume_overlap_with_protein",
    "volume_overlap_with_organic_cofactors",
    "volume_overlap_with_inorganic_cofactors",
    # interreceptor validity
    "protein-ligand_maximum_distance",
]

评价过滤生成分子涉及的一些函数：

def df_split_mol_frags(mol_table: pd.DataFrame) -> pd.DataFrame:
    """Split the molecules in the DataFrame into SDF fragments.


    :param mol_table: Molecule table DataFrame.
    :return: Molecule table DataFrame with fragments.
    """
    new_rows = []
    for row in mol_table.itertuples():
        mols_preds = Chem.SDMolSupplier(str(row.mol_pred), removeHs=False)
        for mol_pred in mols_preds:
            new_row = row._asdict()
            new_row["mol_cond"] = row.mol_cond
            with tempfile.NamedTemporaryFile(suffix=".sdf", delete=False) as temp_pred:
                Chem.MolToMolFile(mol_pred, temp_pred.name)
                new_row["mol_pred"] = temp_pred.name
            new_rows.append(new_row)
    new_df = pd.DataFrame(new_rows)
    new_df["mol_true"] = None
    return new_df




def pb_screen_mols(
    input_molecule_dir: str,
    reference_protein_path: str,
    output_molecule_dir: str,
    bust_results_filepath: str,
):
    """
    Screen a directory of molecule SDF files using the PoseBusters validity metric.
    Only PoseBuster-valid molecule SDF files will be saved to the output directory.


    :param input_molecule_dir: Path to a directory containing SDF files of molecules to screen.
    :param reference_protein_path: Path to a reference protein PDB file.
    :param output_molecule_dir: Path to a directory to which to save PoseBuster-valid SDF files.
    :param bust_results_filepath: Path to which to save the PoseBusters results.
    """
    assert os.path.exists(
        input_molecule_dir
    ), f"Input molecule directory {input_molecule_dir} does not exist."
    assert os.path.exists(
        reference_protein_path
    ), f"Reference protein file {reference_protein_path} does not exist."
    os.makedirs(output_molecule_dir, exist_ok=True)


    # collect inference files


    inference_sdf_results = [
        item
        for item in glob.glob(os.path.join(input_molecule_dir, "*"))
        if item.endswith(".sdf")
    ]
    assert inference_sdf_results, f"No SDF files found in {input_molecule_dir}."
    dataset_protein_files = [
        reference_protein_path for _ in range(len(inference_sdf_results))
    ]
    mol_table = pd.DataFrame(
        {
            "mol_pred": [item for item in inference_sdf_results if item is not None],
            "mol_true": None,
            "mol_cond": [item for item in dataset_protein_files if item is not None],
        }
    )
    mol_frags_table = df_split_mol_frags(mol_table)


    # screen molecules


    buster = PoseBusters(config="dock", top_n=None)
    bust_results = buster.bust_table(mol_frags_table, full_report=False)
    bust_results.to_csv(bust_results_filepath, index=False)
    print(
        f"PoseBusters results for input molecule directory {input_molecule_dir} saved to {bust_results_filepath}."
    )


    # filter invalid molecules


    bust_results_df = pd.read_csv(bust_results_filepath)
    bust_results_df["valid"] = bust_results_df[BUST_TEST_COLUMNS].all(axis=1)
    valid_mol_indices = bust_results_df[bust_results_df["valid"] == 1].index
    valid_mol_paths = mol_frags_table.iloc[valid_mol_indices]["mol_pred"]
    for valid_mol_path in valid_mol_paths:
        shutil.move(
            valid_mol_path,
            os.path.abspath(
                os.path.join(output_molecule_dir, os.path.basename(valid_mol_path))
            ),
        )
    print(f"PoseBuster-valid molecules saved to {output_molecule_dir}.")




if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    parser.add_argument(
        "--input_molecule_dir",
        type=str,
        default=os.path.join(
            "experiments", "4OZ2_1Z6_07_25_2024_12_00_00", "outputs_07_25_2024_12_00_00"
        ),
        help="Path to a directory containing SDF files of molecules to screen.",
    )
    parser.add_argument(
        "--reference_protein_path",
        type=str,
        default=os.path.join(
            "experiments", "4OZ2_1Z6_07_25_2024_12_00_00", "4OZ2_without_1Z6.pdb"
        ),
        help="Path to a reference protein PDB file.",
    )
    parser.add_argument(
        "--output_molecule_dir",
        type=str,
        default=os.path.join(
            "experiments", "4OZ2_1Z6_07_25_2024_12_00_00", "pb_valid_outputs_07_25_2024_12_00_00"
        ),
        help="Path to a directory to which to save PoseBuster-valid SDF files.",
    )
    parser.add_argument(
        "--bust_results_filepath",
        type=str,
        default=os.path.join(
            "experiments",
            "4OZ2_1Z6_07_25_2024_12_00_00",
            "outputs_07_25_2024_12_00_00",
            "bust_results.csv",
        ),
        help="Path to which to save the PoseBusters results.",
    )
    args = parser.parse_args()


    pb_screen_mols(
        args.input_molecule_dir,
        args.reference_protein_path,
        args.output_molecule_dir,
        args.bust_results_filepath,
    )

通过 PoseBusters 评价生成分子，具体命令如下：

python ./experiments/4OZ2_1Z6_07_25_2024_12_00_00/pb_screen_mols.py \
  --input_molecule_dir ./experiments/4OZ2_1Z6_07_25_2024_12_00_00/outputs_07_25_2024_12_00_00 \
  --reference_protein_path ./experiments/4OZ2_1Z6_07_25_2024_12_00_00/4OZ2_without_1Z6.pdb \
  --output_molecule_dir ./experiments/4OZ2_1Z6_07_25_2024_12_00_00/pb_valid_outputs_07_25_2024_12_00_00 \
  --bust_results_filepath ./experiments/4OZ2_1Z6_07_25_2024_12_00_00/outputs_07_25_2024_12_00_00/bust_results.csv

--input_molecule_dir 指定生成分子所在的文件夹；

--reference_protein_path 指定参考配体 1Z6 的结构文件；

--output_molecule_dir 指定评价过程中生成分子产生的缓存结构的文件夹；

--bust_results_filepath 保存每个生成分子在各项指标的通过情况，这里保存为 bust_results.csv 表格。

评估结果保存在 bust_results.csv 表格中，部分结果展示如下，遍历所有生成分子，通过对应测试的对应 True，没有通过的对应 False 。bust_results.csv 部分内容示例：

我们创建 ./experiments/4OZ2_1Z6_07_25_2024_12_00_00/get_PB_Valid_num.ipynb 统计在生成的 84 个分子中，通过全部 19 项测试的分子数目，代码如下：

import pandas as pd

res_df = pd.read_csv('./outputs_07_25_2024_12_00_00/bust_results.csv')

# 统计所有列为 True 的行数
true_rows_count = res_df[res_df.all(axis=1)].shape[0]
print(true_rows_count)

PB中分子有效性的检测部分，生成有效率有85%，好于文章中的 PB 有效率 46% （可能是因为使用参考分子的结果）。文章中 PB 有效率指的是分子是合理的，包括：sanitizable， 全原子连接的，价键正确，键长键角正确，没有内部碰撞，芳香环共平面等。一共有 21 个分子通过 PoseBusters 的全部测试，占比为 25% 。没有通过全部测试的主要原因是生成分子和蛋白/水的最小距离（minimum_distance_to_protein 或者 minimum_distance_to_waters）较近。

3.2.1.3 评估生成分子的类药属性

创建 ./get_mol_metrics.ipynb 来计算生成分子 QED（类药性）、SA（合成可及性）、logP、lipinski 五规则 和 diversity（多样性）等指标。具体代码如下：

from rdkit import Chem
from analysis.metrics import MoleculeProperties
mol_metrics = MoleculeProperties()

generated_mols = Chem.SDMolSupplier('./experiments/4OZ2_1Z6_07_25_2024_12_00_00/outputs_07_25_2024_12_00_00/4OZ2_mol.sdf')
pocket_mols = [generated_mols]
all_qed, all_sa, all_logp, all_lipinski, per_pocket_diversity = mol_metrics.evaluate(pocket_mols)

运行输出：

84 molecules from 1 pockets evaluated.
QED: 0.539 \pm 0.12
SA: 0.662 \pm 0.09
LogP: 0.737 \pm 1.37
Lipinski: 4.952 \pm 0.26
Diversity: 0.834 \pm 0.00

生成的 84 个分子的 QED 、 SA 和 LogP 分别为 0.539、0.662 和 0.737。分子平均满足 Lipinski 规则的比重较高，均值为 4.952（最高为 5），分子多样性为 0.834 。生成分子的整体类药属性表现良好。

3.2.1.4 对接生成分子并展示

接下来，我们使用 qvina 把生成分子对接到蛋白口袋中，使用 MOE 把配体 1Z6 和蛋白结构分别保存到 ./experiments/4OZ2_1Z6_07_25_2024_12_00_00/4oz2_ligand.sdf 和 ./experiments/4OZ2_1Z6_07_25_2024_12_00_00/4oz2_protein.pdb。把蛋白结构文件格式从 pdb 转换为 pdbqt 。命令如下：

conda activate adt
cd ./experiments/4OZ2_1Z6_07_25_2024_12_00_00
prepare_receptor4.py -r 4oz2_protein.pdb

创建 ./experiments/4OZ2_1Z6_07_25_2024_12_00_00/1_qvina_docking.ipynb 处理生成分子结构数据，生成对接的脚本。具体代码如下：

导入所需的依赖库

import os
from rdkit import Chem

创建 ./experiments/4OZ2_1Z6_07_25_2024_12_00_00/pdbqt_mol 文件夹，把每个生成分子保存成单个的 SDF 格式文件。

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

all_mols = Chem.SDMolSupplier('./outputs_07_25_2024_12_00_00/4OZ2_mol.sdf')

for i,mol in enumerate(all_mols):
    SDWriter = Chem.SDWriter(f'./pdbqt_mol/{i}.sdf')
    SDWriter.write(mol)

根据参考配体的坐标确定对接中心的坐标，把生成分子的格式从 SDF 转换成 pdbqt 。把分子对接命令写入到 docking.sh ，对接设置在边长 20 Å 的正方体范围内。

all_sdfs = [f'/workspace/xxxx/projects/GCDM-SBDD/experiments/4OZ2_1Z6_07_25_2024_12_00_00/pdbqt_mol/{i}.sdf' for i in range(len(all_mols))]

ligand = Chem.SDMolSupplier('./4oz2_ligand.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/GCDM-SBDD/experiments/4OZ2_1Z6_07_25_2024_12_00_00/pdbqt_mol/{i}.pdbqt
    """
        os.system(commands)
        f.write(f'''qvina2 \
        --receptor /workspace/xxxx/projects/GCDM-SBDD/experiments/4OZ2_1Z6_07_25_2024_12_00_00/4oz2_protein.pdbqt \
        --ligand /workspace/xxxx/projects/GCDM-SBDD/experiments/4OZ2_1Z6_07_25_2024_12_00_00/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/GCDM-SBDD/experiments/4OZ2_1Z6_07_25_2024_12_00_00/pdbqt_mol/{i}.log\n''')

运行对接的脚本，命令如下：

conda activate adt
cd /workspace/xxxx/projects/GCDM-SBDD/experiments/4OZ2_1Z6_07_25_2024_12_00_00
sh docking.sh > docking.log 

每个对接分子保存其对应的记录文件在./experiments/4OZ2_1Z6_07_25_2024_12_00_00/pdbqt_mol 中。创建 ./experiments/4OZ2_1Z6_07_25_2024_12_00_00/2_get_best_aff.ipynb，根据对接结果，把最好的对接打分作为分子属性，升序排列保存到一个 SDF 文件中。

import os
from rdkit import Chem

all_logs = [i for i in os.listdir('/workspace/xxxx/projects/GCDM-SBDD/experiments/4OZ2_1Z6_07_25_2024_12_00_00/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'./pdbqt_mol/{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))

生成分子添加 qvina_score 属性之后，保存到 ./experiments/4OZ2_1Z6_07_25_2024_12_00_00/Generated_molecules.sdf 。所有生成的分子如下：

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

qvina score 排名前 3 的分子在口袋中的 Pose 如下，蓝色的是参考分子，绿色的是生成分子：

3.2.2 自定义的测试案例

我们使用 3wze 的蛋白作为自定义测试案例。创建 ./experiments/3WZE 文件夹，从 PDB 数据库下载的蛋白结构 3wze.pdb 。配体在口袋中的构象如下：

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

3.2.2.1 分子生成

使用 MOE 把配体结构和蛋白结构分别命名为 3wze_ligand.sdf 和 3wze_protein.pdb，保存到 ./experiments/3WZE 文件夹。

创建案例的输出文件夹 ./experiments/3WZE/outputs。创建命令如下：

mkdir -p ./experiments/3WZE/outputs/

通过 MOE 把配体-蛋白复合物结构 3WZE.pdb 中的参考配体 BAX 删除，保存为 3WZE_without_BAX.pdb。3WZE.pdb 和 3WZE_without_BAX.pdb 两个结构文件放在 ./experiments/3WZE 文件夹中，供后续评价生成和评价分子使用。

基于 3WZE 的蛋白口袋生成一批适应配体 BAX 结合位点的 3D 分子具体命令如下：

python3 generate_ligands.py \
  checkpoints/bindingmoad_ca_cond_gcpnet.ckpt \
  --pdbfile experiments/3WZE/3WZE.pdb \
  --outdir experiments/3WZE/outputs/ \
  --ref_ligand A:1201 \
  --n_samples 100 \
  --num_nodes_lig 32 \
  --sanitize \
  --relax

指定模型为 checkpoints/bindingmoad_ca_cond_gcpnet.ckpt ；

--pdbfile experiments/3WZE/3WZE.pdb 指定下载的蛋白质-配体复合物位置；

--outdir experiments/3WZE/outputs/ 指定生成分子的保存路径；

--ref_ligand A:1201 指定参考配体 BAX 的信息，包括所属的 A 链，index 为 1201 ；

--n_samples 100 指定采样生成 100 个分子；

--num_nodes_lig 32 指定参考配体的节点数目，即非氢重原子数目；

--sanitize 表示确保生成分子的结构有效性，检查分子是否符合化学规则，包括环的闭合、化学价键等；

--relax 表示使用 UFF 力场优化生成分子 200 步。

生成分子保存为 ./experiments/3WZE/outputs/3WZE_mol.sdf，设置生成 100 个分子，实际生成 55 个分子。batch_size 和参考配体 BAX 的非氢重原子数目一样，为 32 。生成过程显存占用约 15 GB，用时约 38 分钟。

3.2.2.2 通过 PoseBusters 评估生成分子

使用《3.2 内置的案例测试》创建的 ./experiments/4OZ2_1Z6_07_25_2024_12_00_00/pb_screen_mols.py 脚本，通过 PoseBusters 筛选生成分子。具体命令如下：

python ./experiments/4OZ2_1Z6_07_25_2024_12_00_00/pb_screen_mols.py \
  --input_molecule_dir ./experiments/3WZE/outputs \
  --reference_protein_path ./experiments/3WZE/3WZE_without_BAX.pdb \
  --output_molecule_dir ./experiments/3WZE/pb_valid_outputs \
  --bust_results_filepath ./experiments/3WZE/outputs/bust_results.csv

--input_molecule_dir 指定生成分子所在的文件夹。--reference_protein_path 指定不包含参考配体 BAX 的结构文件。--output_molecule_dir 指定评价过程中生成分子产生的缓存结构的文件夹。--bust_results_filepath 保存每个生成分子在各项指标的通过情况，这里保存为 bust_results.csv 表格。

评估结果保存在 bust_results.csv 表格中，部分结果展示如下，遍历所有生成分子，通过对应测试的对应 True，没有通过的对应 False 。

mol_pred_loaded,mol_cond_loaded,sanitization,all_atoms_connected,bond_lengths,bond_angles,internal_steric_clash,aromatic_ring_flatness,double_bond_flatness,internal_energy,protein-ligand_maximum_distance,minimum_distance_to_protein,minimum_distance_to_organic_cofactors,minimum_distance_to_inorganic_cofactors,minimum_distance_to_waters,volume_overlap_with_protein,volume_overlap_with_organic_cofactors,volume_overlap_with_inorganic_cofactors,volume_overlap_with_waters
True,True,True,True,True,True,True,True,True,True,True,False,True,True,False,True,True,True,True
True,True,True,True,True,True,True,True,True,True,True,False,True,True,True,False,True,True,True
True,True,True,True,True,True,True,True,True,True,True,False,True,True,True,False,True,True,True
True,True,True,True,True,True,True,True,True,True,True,False,True,True,True,False,True,True,True

我们创建 ./experiments/3WZE/get_PB_Valid_num.ipynb 统计在生成的 55 个分子中，通过全部 19 项测试的分子数目，代码如下：

import pandas as pd

res_df = pd.read_csv('./outputs/bust_results.csv')

# 统计所有列为 True 的行数
true_rows_count = res_df[res_df.all(axis=1)].shape[0]
print(true_rows_count)

一共有 2 个分子通过 PoseBusters 的全部测试，占比仅为 3.6% 。生成分子不能通过测试的两个主要原因是：1. 有些生成分子和蛋白之间有空间冲突；2. 有些生成分子和蛋白之间没有空间冲突，但距离过近。

3.2.2.3 评估生成分子的类药属性

在 ./get_mol_metrics.ipynb 中计算生成分子 QED（类药性）、SA（合成可及性）、logP、lipinski 五规则 和 diversity（多样性）等指标。

具体代码如下：

from rdkit import Chem
from analysis.metrics import MoleculeProperties
mol_metrics = MoleculeProperties()

generated_mols = Chem.SDMolSupplier('./experiments/3WZE/outputs/3WZE_mol.sdf')
pocket_mols = [generated_mols]
all_qed, all_sa, all_logp, all_lipinski, per_pocket_diversity = mol_metrics.evaluate(pocket_mols)

运行输出：

55 molecules from 1 pockets evaluated.
QED: 0.442 \pm 0.18
SA: 0.539 \pm 0.06
LogP: 1.904 \pm 1.91
Lipinski: 4.691 \pm 0.57
Diversity: 0.723 \pm 0.00

生成的 55 个分子的 QED 、 SA 和 LogP 分别为 0.442、0.539 和 1.904。分子平均满足 Lipinski 规则的比重较高，均值为 4.691（最高为 5），分子多样性为 0.723 。生成分子的类药属性表现一般。

3.2.2.4 对接生成分子并展示

接下来，我们使用 qvina 把生成分子对接到蛋白口袋中，使用 MOE 把配体 BAX 和蛋白结构分别保存到 ./experiments/3WZE/3wze_ligand.sdf 和 ./experiments/3WZE/3wze_protein.pdb。把蛋白结构文件格式从 pdb 转换为 pdbqt 。命令如下：

conda activate adt
cd ./experiments/3WZE
prepare_receptor4.py -r 3wze_protein.pdb

创建 ./experiments/3WZE/1_qvina_docking.ipynb 处理生成分子结构数据，生成对接的脚本。具体代码如下：

导入所需的依赖库

import os
from rdkit import Chem

创建 ./experiments/3WZE/pdbqt_mol 文件夹，每个生成分子保存成单个的 SDF 格式文件。

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

all_mols = Chem.SDMolSupplier('./outputs/3WZE_mol.sdf')

for i,mol in enumerate(all_mols):
    SDWriter = Chem.SDWriter(f'./pdbqt_mol/{i}.sdf')
    SDWriter.write(mol)

根据参考配体的坐标确定对接中心的坐标，把生成分子的格式从 SDF 转换成 pdbqt 。把分子对接命令写入到 docking.sh ，对接设置在边长 20 Å 的正方体范围内。

all_sdfs = [f'/workspace/xxxx/projects/GCDM-SBDD/experiments/3WZE/pdbqt_mol/{i}.sdf' for i in range(len(all_mols))]

ligand = Chem.SDMolSupplier('./3wze_ligand.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/GCDM-SBDD/experiments/3WZE/pdbqt_mol/{i}.pdbqt
    """
        os.system(commands)
        f.write(f'''qvina2 \
        --receptor /workspace/xxxx/projects/GCDM-SBDD/experiments/3WZE/3wze_protein.pdbqt \
        --ligand /workspace/xxxx/projects/GCDM-SBDD/experiments/3WZE/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/GCDM-SBDD/experiments/3WZE/pdbqt_mol/{i}.log\n''')

运行对接的脚本，命令如下：

conda activate adt
cd /workspace/xxxx/projects/GCDM-SBDD/experiments/3WZE
sh docking.sh > docking.log 

每个对接分子保存其对应的记录文件在 ./experiments/3WZE/pdbqt_mol 中。创建 ./experiments/3WZE/2_get_best_aff.ipynb，根据对接结果，把最好的对接打分作为分子属性，升序排列保存到一个 SDF 文件中。

import os
from rdkit import Chem

all_logs = [i for i in os.listdir('/workspace/xxxx/projects/GCDM-SBDD/experiments/3WZE/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'./pdbqt_mol/{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))

生成分子添加 qvina_score 属性之后，保存到 ./experiments/3WZE/Generated_molecules.sdf 。所有生成的分子如下：

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

qvina score 排名前 3 的分子在口袋中的 Pose 如下，绿色的是参考分子，紫红色的是生成分子。可以看出生成分子和蛋白存在空间冲突问题：

四、总结

GCDM 的主要创新点和意义在于：

(1) 几何完备性：GCDM 首次将几何完备的图神经网络与扩散模型相结合，实现了对 3D 分子结构的精确建模和生成。这一创新使得模型能够捕获分子的完整几何信息，包括手性特性。

(2) 大分子生成能力: 通过引入几何完备的表示学习，GCDM 成功克服了现有方法在生成大型3D 分子时的局限性。这一突破为复杂药物分子的设计和优化提供了新的可能。

(3) 多任务适用性: GCDM 在无条件生成、条件生成、分子优化和蛋白质靶向分子设计等多个任务上均表现出色，展现出广泛的应用潜力。

(4) 计算效率: 尽管 GCDM 引入了更复杂的几何表示，但其计算效率仍然可观。例如，生成 250 个新的大分子仅需约 15 分钟（GPU 为 24GB NVIDIA A10）。

(5) 可解释性: GCDM 生成的分子结构具有良好的物理和化学合理性，这为进一步理解分子设计过程提供了洞察。

虽然作者提出了一个很好的模型，但是实际测试下来，在有口袋条件下，生成的分子，根据视频结果来看，生成分子的质量比较糟糕，口袋的限制也几乎没有。可能是作者提供了一个未充分训练的 checkpoint，具体还需要从头训练模型才知道。另一方面，这些生成的分子质量那么差，类药性，合格率都很高，这说明，这些计算指标并不真实。