DeepICL 是一个基于相互作用感知的 3D 分子生成模型，能够在目标结合口袋内进行相互作用引导的小分子设计。DeepICL 通过利用蛋白质-配体相互作用的普遍模式作为先验知识，在有限的实验数据下也能实现高度的泛化能力。

一、背景介绍

DeepICL 来源于韩国科学技术院化学系和人工智能研究所的 Woo Youn Kim 教授为通讯作者的文章：《3D molecular generative framework for interaction-guided drug design》。该文章于 2024 年 3 月 27 日发表在 《Nature Communications》上。文章链接：​​​​​​​​​​​​​​3D molecular generative framework for interaction-guided drug design | Nature Communications。

当前基于结构的分子生成模型由于受限于缺少数据，在泛化方面面临挑战，对没有学习过的目标蛋白质往往产生不利的相互作用。为了解决这一问题，作者提出了 DeepICL ，学习蛋白质-配体相互作用的普遍模式并将其作为先验知识，可以在有限的实验数据下也能实现高度的泛化能力。通过分析生成的配体在未见过的目标中的结合位姿稳定性、亲和力、几何模式、多样性和新颖性，作者对 DeepICL 模型性能进行了全面评估。此外，潜在突变体选择性抑制剂的有效设计也证明了DeepICL 在基于结构的药物设计中的适用性。

二、模型介绍

基于结构的分子生成方法由于缺少活性数据，泛化能力受到限制，生成的分子不够新颖，还往往和目标蛋白之间产生不利的相互作用。整合足够的先验知识有助于开发具有泛化能力的深度学习模型。

为了解决这个问题，作者提出了 DeepICL （Deep Interaction-aware Conditional Ligand generative model），利用蛋白质-配体相互作用的普遍性质来引导基于结构的药物设计。

DeepICL 模型使用配体应该结合的子口袋中的局部相互作用条件，而不是使用整个相互作用信息，以防止对特定口袋或配体结构产生不良偏向。

2.1 模型框架

模型框架包括两个主要阶段：（1）相互作用感知条件设置；（2）相互作用感知的三维分子生成。

框架的第一阶段旨在通过研究给定结合位点的蛋白质原子来设定相互作用条件。作者使用了四种类型的蛋白质-配体相互作用：氢键、盐桥、疏水相互作用和π–π堆积。作者主要通过两种策略确定口袋原子的相互作用类别，如下图 a 所示。第一种策略，基于口袋原子定义相互作用，使用在在生成阶段，受体如何与配体相互作用的信息并不是预知的。因此，作者预定义了相互作用类别的标准（SMART 模式），以便通过分析来指定每个蛋白质原子的相互作用条件。由于没有使用任何参考配体来设置条件，作者称之为无参考相互作用条件。第二种，使用PLIP 提取相互作用，使用在在训练阶段期间，有蛋白质-配体复合物的真实结构，作者使用蛋白质-配体相互作用分析器（PLIP）从这些参考结构中提取相互作用条件。除了这两种策略外，还可以根据对目标系统的理解手动指定所需的相互作用条件。

关于相互作用，作者将相互作用条件定义为一组蛋白质原子的附加相互作用类型one-hot 向量，它表示原子是否可以参与特定相互作用及其在相互作用中的作用。蛋白质原子分为七类之一: 阴离子、阳离子、氢键供体和受体、芳香族、疏水和非相互作用原子。与将整个相互作用信息表示为单个相互作用指纹相比，这一策略旨在局部建立相互作用条件。

在第二阶段，基于口袋的三维上下文和第一阶段的相互作用条件，逐步生成配体中的原子。如下图 b 所示，配体原子一个接一个地添加，且每一步的兴趣原子发生变化。这个兴趣原子指示下一原子添加的位置。因此，作者只考虑兴趣原子C_{t}的周围环境，以便重新定义局部相互作用条件I_{t}并输入到当前生成步骤 t 中。作者针对目标蛋白设计了两个分子生成任务的应用——配体优化（Ligand elaboration）和从头配体设计（de novo）。前者优化配体结构以提高其药效，后者从头设计配体，可以生成适应口袋的多样化分子。对于配体优化任务，给出了配体核心结构的结合 pose 并用作初始状态。在从头配体设计任务中，可以手动选择口袋中的一个点作为起始点。

DeepICL 模型是由编码器和解码器组成的变分自编码器 VAE 框架。如下图所示，编码器 (Encode) 模块将给定的蛋白质-配体复合体的结构嵌入到一个遵循标准正态分布的潜在向量 z 中。然后，解码器（Decode）模块从潜在向量 z 中以逐原子的方式顺序生成配体结构。相互作用条件 I 被整合到潜在向量 z 中，以放置一个适合的配体原子以与目标形成期望的相互作用。编码器和解码器模块共享相同的嵌入层，这些嵌入层由多层E(3)-不变的相互作用网络组成，用于在蛋白质和配体之间传播消息。

DeepICL 使用两个额外的虚拟原子，它们仅持有位置信息，即质心和关注原子，来辅助配体设计过程。原始配体的质心大致确定了待生成配体的全局位置。关注原子限制了下一个配体原子将被放置的3D空间；在每一步中，仅考虑其邻近的蛋白质原子来预测下一个原子的类型及其位置。因此，DeepICL能够学习局部口袋环境与配体结构偏好之间的关系，以满足给定的相互作用条件，利用DeepICL在任何蛋白质的配体设计任务中的鲁棒性。在训练和采样过程中，上述两个虚拟原子被作为单独的配体原子对待，然后在最终确定配体结构时被移除。

2.2 训练和测试数据

在这项工作中，作者仅使用了 PDBbind 2020 中的 general set 里的实验晶体结构，这些结合结构是通过 X 射线晶体学鉴定的。作者根据目标序列的相似性对晶体结构数据进行了划分，以确保训练和测试数据之间的任何数据对的相似性都不超过 0.6，这些相似性是通过 CD-HIT 软件计算和聚类得到。数据处理后，作者使用了 11284 个结构来训练模型，2109 个结构用于验证。作者筛选出了剩余的数据，保留了 109 个满足以下三个条件的测试复合体：(1) 配体的Tanimoto 相似性与训练集中所有配体的相似性都小于 0.6；(2) 每条数据对应不同的蛋白质家族；(3) 蛋白质的重原子数少于 3000。为了方便起见，作者从中随机选择了100个测试复合体。

2.3 模型性能

2.3.1 基于相互作用条件的配体设计效果

作者首先展示了相互作用感知条件在设计特定相互作用模式上的效果。参考蛋白质-配体化合物相互作用模式，来评估生成分子在相互作用模式方面的表现。作者选择了三种蛋白质，分别是骨形态发生蛋白1（BMP1）、成纤维细胞生长因子1（FGF1）和二氢叶酸还原酶（DHFR）。从它们原始复合物的参考配体中，提取了初始核心结构。核心结构是参考配体去除了链和官能团，保留了由单环或双环组成的最小结构，如下图 a 所示。从下图 b 可以看出，设计的新配体和参考配体具有很高的相互作用相似性。

下图 c 和 d 分别展示了原始配体和设计的配体（具有最高相互作用相似性）的相互作用图。对于第一个案例，模型成功设计了氢键、π–π 堆积和盐桥，符合给定的相互作用条件但具有不同的结构特征。值得注意的是，下图 d-1 显示模型能够生成噻吩环代替原始的苯环，以与酪氨酸 68（TYR68）形成 π–π 堆积。这表明模型学习到了芳香族基团所需的特征以形成 π–π 堆积。虽然模型在疏水的苯丙氨酸 157（PHE157）附近添加了脂肪族碳，但距离略大于被概况为疏水相互作用的阈值。对于其他两个案例，DeepICL也成功设计了显示与原始配体高度相似的相互作用模式的配体，同时生成了与原始配体不同的结构特征。

进一步证明了相互作用感知条件策略在 DeepICL 的效果，作者设计了消融实验，屏蔽相互作用条件以排除推断过程中的参考相互作用模式信息。下图中比较了使用参考条件或屏蔽条件生成的两个配体集合之间的相互作用相似性分布 (相互作用指纹的 cosine similarity)。很明显，每个案例中使用相互作用参考条件扩展的配体具有更高的相互作用相似性，从而证明了DeepICL 框架是高度可控的，通过扩展现有配体可以提供具有期望相互作用模式的配体。

2.3.2 展示框架的通用性

在后续的评估中，为了考察相互作用感知生成框架 DeepICL 在基于结构的药物设计中的通用性，作者设计了一个只在结合结构上训练的模型（基线模型，即图中的 Without information），没有任何关于蛋白质-配体相互作用的明确信息，更具体地说，没有使用用于蛋白质原子的附加相互作用条件向量。基线模型可能会基于原子出现频率学习一些与非共价相互作用相关的信息，但由于训练集中蛋白质-配体对的数量有限，不太可能在典型的非共价相互作用模式上被泛化。因此，基线模型在确定新添加原子的类型和位置时，不可避免地依赖于蛋白质-配体结合几何的统计分布。在此，作者将有无相互作用信息生成的配体集合分别命名为交互作用条件模型和基线模型。

2.3.3 结合位姿稳定性分析

为评估设计的配体的结合稳定性，作者使用基线模型也设计了新配体，随机选择了和参考配体重原子数目相同的 10 个设计配体进行 MD 模拟并统计了 RMSD 变化。从下图 a 绘制了 10 个设计配体的平均 RMSD 值以及原始配体在模拟过程中的 RMSD 值。值得注意的是，基于相互作用信息的配体显示出的 RMSD 值与原始配体相当低。这是模型泛化能力的有力证据，即模型能够生成与参考配体一样有效地稳定结合于未见过靶点的配体。尽管模型未在显式指示结合稳定性的配体 MD 轨迹数据上进行训练，但它成功地建立了稳定的相互作用。此外，未使用相互作用信息模型生成的配体在所有三个案例中显示了更高的 RMSD 值，这暗示了缺乏稳定其结合位姿所需的有利相互作用。在 BMP1 和 DHFR 的案例中，两组之间显示出了明显的差异，其中未使用相互作用信息的配体与其初始结合位姿的偏差显著增大。同时，FGF1 的两个集合之间差异相对较小，两者相比于原始配体都显示出可接受的稳定性。这个结果表明，即使是基线模型，也可以在一定程度上通过在训练阶段学习蛋白质-配体相互作用来生成稳定的配体。

2.3.4 从头设计配体的结合亲和力分析

为了分析从头生成配体的结合亲和力，作者使用 100 个测试蛋白分别生成 100 个配体，然后通过 SMINA 评估结合亲和力。下图 b 展示了配体结合亲和力分数的统计数据。基线模型生成的配体的平均结合亲和力值为 -6.52 kcal/mol，高于使用相互作用信息生成的配体（-7.67 kcal/mol）。这一结果表明结合相互作用信息有助于设计具有更强结合亲和力的分子。作者的的策略通过在没有实验亲和力数据训练的情况下实现相当高的结合亲和力，帮助提高了生成模型的泛化能力。

作者还分析了每种类型的蛋白质-配体相互作用的数量，以解释相互作用信息在实现高结合亲和力中的作用。由于没有具体的期望相互作用数量的值，作者比较了在相同测试口袋下有无相互作用信息生成的复合物。如下图 c 所示。对于基线模型，疏水相互作用和氢键的数量远少于使用相互作用信息训练的模型，而盐桥和 π–π 堆积相当。没有相互作用信息生成的每个分子的氢键数量仅为有信息的约一半，这一差异远大于其他类型的相互作用。氢键是供体和受体之间的方向性相互作用，当仅依赖于有限数量的结构数据的分布而没有蛋白质-配体相互作用的先验知识时，生成这种相互作用类型更加具有挑战性。（换句话说，如果不使用相互作用作为条件训练模型，那么模型生成的分子产生的H见，和疏水相互作用比较少）

与基线模型相比，DeepICL 通过在部分确定性方式下采样下一个原子（由作为条件给定的蛋白质-配体相互作用的先验知识指导）可以生成更多的氢键原子对。在生成氢键受偏爱情况下，使用相互作用信息的模型添加 N、O 和 F 的比例高于基线模型。作者将使用相互作用信息生成的配体更高的结合亲和力归因于氢键形成的更高成功率。

2.3.5 生成相互作用的几何分析

大多数当前基于结构的药物设计的深度生成建模方法仅关注设计的配体的分子内几何形态，而忽略了分子间几何的分析。因此，作者展示了在采样的复合物中蛋白质-配体相互作用的几何形态。对于结合亲和力分析中生成的配体，作者测量了每种非共价相互作用类型的距离，而没有进行进一步的结构优化。下图 d 和 e 展示了疏水相互作用和氢键的几何分布。下图 d 显示了疏水相互作用距离的密度分布，这种类型在 PDB 结构中最为常见，DeepICL 有效捕捉了随着距离减小密度衰减的观察趋势。下图 e 展示了氢键距离的密度分布。氢键的距离约为 3.0Å ，生成数据的分布也显示出一个接近 3.0 Å 的峰值。在这两种类型中，无论是否使用相互作用信息，分布都没有变化。未经过显式相互作用信息训练的基线模型也能够捕捉这些相互作用类型的空间分布。这表明关于化学相互作用的知识有助于模型在特定结合点生成正确的相互作用类型，而不是形成更合理的几何模式，从而提高了有利相互作用的形成率。

2.3.6 设计配体的化学多样性和新颖性

实现高化学多样性和新颖性是基于结构的药物设计中的另一个重要目标，这可以在核心结构（或称为骨架）层面进行评估。作者评估了 100 个测试蛋白生成的 10000 个配体的骨架唯一性。在9930 个有效分子中，去除重复项后得到了 5669 个唯一骨架，骨架多样性为 57.1%。为了比较，作者还评估了训练数据的骨架多样性。在 10752 个分子中，拥有 5783 个唯一骨架，骨架多样性为 53.8% 。值得注意的是，尽管使用了从参考配体中提取的特定相互作用条件，DeepICL 仍实现了略高于训练数据的多样性。下图 a 展示了生成配体的前十个骨架频率情况。模型最常生成苯环，它也是训练数据中最常见的骨架。模型还生成了训练数据中不太常见的骨架分子；例如，二苯基甲烷在生成的配体中排名第六，而在训练配体中排名第十七。萘则表现出更大的差异：在生成的配体中排名第九，而在训练数据中排名第三百五十二。因此，模型并不完全遵循训练数据中观察到的结构优先级。

接着，作者评估了设计配体的结构新颖性。与训练数据相比，在 5669 个唯一骨架中有 5467 个骨架是新颖的，达到了 96.4% 的新颖性。这表明框架能够提供新颖的结构，而不是重复从训练数据中学习到的核心结构。下图 b 展示了一些新颖生成骨架的示例。

2.3.7 针对特定位点的相互作用调控以选择性地控制配体设计

DeepICL 框架的一个关键优势是能够基于已有知识建立相互作用条件，从而设计具有特定功能的配体。在此，作者选择了一个在特定位置形成选择性相互作用至关重要的实际问题：设计一种能够选择性结合突变型表皮生长因子受体（EGFR）而不影响野生型 EGFR 的配体。下图 a 展示了设计配体在野生型和突变型 EGFR 上的结合亲和力分数。位于实线对角线以上的点在突变口袋上的得分低于野生型。较低的分数表明更强的结合，这种趋势清楚地表明生成的配体主要能更强地结合突变型 EGFR 。由于能量减少 1.36 kcal/mol 在理论上相当于抑制浓度降低 10 倍，作者设置 2.72 kcal/mol 的差异或抑制浓度降低 100 倍作为鉴别配体选择性的标准。结果，作者获得了233 个选择性配体，对应于下图 a 中虚线以上的红点。作者在选择性配体中通过目视检查选择了一个设计良好的配体，并在下图 b 中进行了可视化。该配体与 MET790 的侧链形成了疏水相互作用，同时与 ARG858 的主链形成了氢键，成功地通过使用特定位点调控方法识别出一个与突变残基相互作用的有前景的分子，而无需额外数据进行训练。

三、DeepICL 评测

3.1 安装环境

复制代码项目：

git clone https://github.com/ACE-KAIST/DeepICL.git

根据项目提供的安装脚本（install_packages.sh）创建 DeepICL 环境

chmod +x install_packages.sh
bash install_packages.sh

脚本中包含安装的依赖包以及版本信息，具体内容如下：

#!/bin/bash


# 1. Create conda environment
conda create -n DeepICL python=3.7.13 -y
conda activate DeepICL


# 2. install packages dependent on mkl
conda install scipy=1.6.2 numpy=1.21.2 pandas=1.3.4 scikit-learn=1.0.2 seaborn=0.11.0 -y


# 3. install pytorch and torch_geometric
conda install pytorch=1.9.0 cudatoolkit=10.2 -c pytorch -y
conda install pyg=2.0.3 -c pyg -c conda-forge -y


# 4. install others
conda install -c rdkit rdkit=2022.09.1 -y
conda install -c conda-forge biopython=1.77 openbabel=3.1.1 -y
conda install -c conda-forge plip=2.2.2 -y

3.2 数据预处理

模型使用的数据集是 PDB V2020 的 general set ，链接是 Welcome to PDBbind-CN database 。点击下面红色方框位置下载原始数据集。

下载好的 PDBbind_v2020_other_PL.tar.gz 放在项目文件夹 ./ 目录中，解压命令如下：

tar -xzvf PDBbind_v2020_other_PL.tar.gz

解压后的文件夹为 v2020-other-PL，方便起见重命名为 v2020 。此时，完整的 DeepICL 项目的文件结构如下：

.
|-- LICENSE
|-- PDBbind_v2020_other_PL.tar.gz
|-- README.md
|-- crossdocked
|-- data
|-- example
|-- image
|-- install_packages.sh
|-- processed_data
|-- script
|-- temp
|-- test
`-- v2020


9 directories, 4 files

对原始数据集进行处理，命令如下：

python ./data/preprocessing.py \
  ./v2020 \
  ./processed_data \
  128

使用 ./data/preprocessing.py 处理原始蛋白质-配体结构数据，创建 ./processed_data 文件夹以保存处理好的数据。第一个参数 ./v2020 指定下载的 PDB v2000 的 general set 的原始结构数据。第二个参数 ./processed_data 指定处理好的数据的保存路径。第三个参数 128 表示数据处理使用的 CPU 核数，我们使用的机器是 128 个核心的，这里根据机器情况设置。

运行输出示例：

save_dir: ./processed_data
NUM DATA: 14127
Traceback (most recent call last):
  File "/workspace/xxxx/projects/DeepICL/data/processor.py", line 567, in _processor
    assert ligand_mol is not None  # , "ligand mol is None"
AssertionError
 ... 
 
6om4
Exceed max ligand atom num
Traceback (most recent call last):
  File "/workspace/xxxx/projects/DeepICL/data/processor.py", line 579, in _processor
    assert self._filter_ligand(ligand_mol), ligand_fn.split('/')[-1]
AssertionError: 6pek_ligand.sdf


6pek
NUM PROCESSED DATA: 9705
PROCESSING TIME: 233.9 (s)

原始的结构数据一共是 14127 个蛋白，作者把配体原子数目超过 50 个的过滤掉，有些配体无法正常读取。数据处理完用时大约 4 分钟，得到 9705 个处理好的文件，保存在 ./processed_data 中，每个正确处理好的蛋白信息保存为一个二进制的文件, 例如: ./processed_data/1a0q，名字为对应的 PDB ID 。处理过程的记录文件在 ./process.log 

以处理好的 ./processed_data/1a0q 为例，介绍保存的信息，参考 ./check_processed_data.ipynb。

import pickle


# 文件路径
file_path = './processed_data/1a0q'


# 读取数据
with open(file_path, "rb") as r:
    data = pickle.load(r)

data 处理后的文件（例如：./processed_data/1a0q）保存到一个包含 13 个元素的元组中，内容包括 ligand_type , ligand_coord , pocket_type , pocket_coord , ligand_adj , None , None , ligand_n , 0 , ligand_smi , None , pocket_cond , center_of_mass 。

主要内容为：配体和口袋的非氢原子类型和坐标（ligand_type , ligand_coord , pocket_type , pocket_coord ），配体的邻接矩阵（ligand_adj），原子数（ligand_n），smiles（ligand_smi），相互作用信息作为口袋的条件（pocket_cond），配体质心（center_of_mass） 。下面查看主要的内容。

data[0] 为配体的原子类型（ligand_type），ATOM_TYPES = ["C", "N", "O", "F", "P", "S", "Cl", "Br"] 和一个未知项，9 维的 one-hot 向量表示一个原子的类型，ligand 包含 23 个非氢原子，矩阵形状为 （23，9）。具体内容如下：

array([[1., 0., 0., 0., 0., 0., 0., 0., 0.],
       [1., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 1., 0., 0., 0., 0.],
       [0., 1., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 1., 0., 0., 0., 0., 0., 0.],
       [0., 0., 1., 0., 0., 0., 0., 0., 0.],
       [0., 0., 1., 0., 0., 0., 0., 0., 0.],
       [1., 0., 0., 0., 0., 0., 0., 0., 0.],
       [1., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 1., 0., 0., 0., 0., 0., 0.],
       [1., 0., 0., 0., 0., 0., 0., 0., 0.],
       [1., 0., 0., 0., 0., 0., 0., 0., 0.],
       [1., 0., 0., 0., 0., 0., 0., 0., 0.],
       [1., 0., 0., 0., 0., 0., 0., 0., 0.],
       [1., 0., 0., 0., 0., 0., 0., 0., 0.],
       [1., 0., 0., 0., 0., 0., 0., 0., 0.],
       [1., 0., 0., 0., 0., 0., 0., 0., 0.],
       [1., 0., 0., 0., 0., 0., 0., 0., 0.],
       [1., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 1., 0., 0., 0., 0., 0., 0.],
       [1., 0., 0., 0., 0., 0., 0., 0., 0.],
       [1., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 1., 0., 0., 0., 0., 0., 0.]])

data[1] 对应 ligand 的原子坐标（ligand_coord），矩阵形状为 （23，3）。内容如下：

array([[-0.14473832,  0.02331776,  0.18624299],
       [-0.99773832, -0.83468224, -0.74875701],
       ... 
       [ 2.93426168,  1.03331776, -5.47775701],
       [-3.48673832, -1.54268224,  5.09224299]])

data[2] 对应蛋白口袋的原子特征（pocket_type），用一个 51 维的特征向量表示一个原子，包括原子类型(9)、与其他原子的成键数目(6)、原子杂化类型(7)，电荷数目(7)，是否属于芳香结构(1)，氨基酸类型(20)等信息。包含 162 个非氢原子，矩阵形状为 （162，51）。数据示意如下：

array([[0., 1., 0., ..., 0., 0., 0.],
       [1., 0., 0., ..., 0., 0., 0.],
       [1., 0., 0., ..., 0., 0., 0.],
       ...,
       [1., 0., 0., ..., 0., 0., 0.],
       [1., 0., 0., ..., 0., 0., 0.],
       [1., 0., 0., ..., 0., 0., 0.]])

data[3] 对应蛋白口袋的原子坐标（pocket_coord），矩阵形状为 （162，3）。数据示意如下：

array([[-4.17738318e-01, -8.11868224e+00, -2.58775701e+00],
       [-1.42973832e+00, -7.51268224e+00, -3.42775701e+00],
       [-1.20773832e+00, -7.82368224e+00, -4.90775701e+00],
       ...
       [ 2.26261682e-01,  2.00331776e+00, -8.69775701e+00],
       [ 2.15826168e+00,  3.24331776e+00, -9.47275701e+00],
       [ 1.59926168e+00,  2.08031776e+00, -8.93875701e+00]])

data[4] 对应配体的邻接矩阵（ligand_adj），矩阵形状为 （23，23）。数据示意如下：

array([[0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0],
       [1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0],
       ...
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0,
        0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
        0]], dtype=int32)

data[7] 和 data[9] 分别对应配体的原子数目和 SMILES，即 '23' 和 'CCCC[C@H](NC(=O)CCC(=O)O)[P@@](=O)(O)Oc1ccccc1'

data[10] 是小分子 rdkit.Chem.Scaffolds.MurckoScaffold 提取出来的骨架smiles。

data[11] 表示蛋白-配体的相互作用信息，作为口袋的条件（pocket_cond），是通过 plip 库计算得到的。每个原子的相互作用用一个 7 维的 one-hot 向量表示，INTERACTION_TYPES = ["pipi", "anion", "cation", "hbd", "hba", "hydro"]，包括 π-π 相互作用、阴离子相互作用、阳离子相互作用、氢键供体、氢键受体和疏水相互作用等六种相互作用。六种相互作用和一个未知项共同表示一个原子的相互作用信息，对应 7 维的 one-hot 向量表示。相互作用矩阵的形状为（162，7）;162 对应的是口袋的原子数，数据示意如下：

array([[0., 0., 0., ..., 0., 0., 1.],
       [0., 0., 0., ..., 0., 0., 1.],
       [0., 0., 0., ..., 0., 0., 1.],
       ...,
       [0., 0., 0., ..., 0., 0., 1.],
       [0., 0., 0., ..., 0., 0., 1.],
       [0., 0., 0., ..., 0., 0., 1.]])

（注：相互作用 one-hot 编码是添加在口袋的原子上的，与 atom type 类似）

data[12] 对应内容为配体的质心坐标（center_of_mass），即 [13.44973832, 20.67668224, 59.21875701] 。

数据处理完后，在训练之前，需要更新生成训练集和验证集的index文件（train_keys.pkl 和valid_keys.pkl）。使用到我们自己的脚本./sripts/generate_new_key.py。运行方法：

python script/generate_new_key.py \
  ./processed_data/ \
  0.1 \
  ./data/new_keys/

参数中，指定处理好的数据保存在 ./processed_data/；valid的数据比例是0.1; 生成的train和valid pkl文件保存在 ./data/new_keys/。

3.3 训练模型

由于作者没有提供训练好的 checkpoint，因此，在分子生成前需要先训练模型。(注：模型训练部分的代码进行了修改，主要是解决了参数等报错)

使用上述准备好的数据（./processed_data/）以及 index文件（./data/new_keys/）,训练 DeepICL 的命令如下：

python -u ./script/train.py \
  --world_size 1 \
  --save_dir ./save_model/ \
  --data_dir ./processed_data/ \
  --key_dir ./data/new_keys/ \
  --lr 1e-3 \
  --num_epochs 1001 \
  --save_every 1 \
  --k 8 \
  --lr_decay 0.8 \
  --lr_tolerance 4 \
  --lr_min 1e-6 \
  --num_workers 13 \
  --conditional

训练过程输出：

############################ Finished Loading Model ############################
EPOCH || TRA_VAE | TRA_SSL | TRA_TYPE | TRA_DIST | TRA_TOT || VAL_VAE | VAL_SSL | VAL_TYPE | VAL_DIST | VAL_TOT || TIME/EPOCH | LR | BEST_EPOCH
[W reducer.cpp:1158] Warning: find_unused_parameters=True was specified in DDP constructor, but did not find any unused parameters in the forward pass. This flag results in an extra traversal of the autograd graph every iteration,  which can adversely affect performance. If your model indeed never has any unused parameters in the forward pass, consider turning this flag off. Note that this warning may be a false positive if your model has flow control causing later iterations to have unused parameters. (function operator())
0 || 0.000 | 0.000 | 1.400 | 2.821 | 4.221 || 0.000 | 0.000 | 1.125 | 2.266 | 3.391 || 2040.89 | 0.0010 | 0*
1 || 0.027 | 0.000 | 1.038 | 2.126 | 3.190 || 0.001 | 0.000 | 0.979 | 2.026 | 3.006 || 1962.47 | 0.0010 | 1*

注：因为作者使用 pytorch 的数据并行，即使我们只有一张显卡，但是还是占用了大量的CPU资源。此外，模型在训练时，CPU和 GPU的占用，但是CPU占用率比较高近乎100%，GPU的占用率较低只有30%左右。训练了4个小时，只有8个epoches，每个epoch需要33分钟。。

训练过程中，每一个epoch保存一个checkpoint，保存在路径：./save_model。

但是训练到一般就出现了 nan，如下：

13 || 0.005 | 0.000 | 1.062 | 2.091 | 3.158 || 0.001 | 0.000 | 0.915 | 1.795 | 2.711 || 1886.38 | 0.0005 | 7
14 || nan | 0.000 | nan | nan | nan || nan | 0.000 | nan | nan | nan || 1884.46 | 0.0004 | 7
15 || nan | 0.000 | nan | nan | nan || nan | 0.000 | nan | nan | nan || 1881.60 | 0.0003 | 7
16 || nan | 0.000 | nan | nan | nan || nan | 0.000 | nan | nan | nan || 1879.88 | 0.0003 | 7
17 || nan | 0.000 | nan | nan | nan || nan | 0.000 | nan | nan | nan || 1883.93 | 0.0002 | 7
18 || nan | 0.000 | nan | nan | nan || nan | 0.000 | nan | nan | nan || 1885.69 | 0.0002 | 7
19 || nan | 0.000 | nan | nan | nan || nan | 0.000 | nan | nan | nan || 1884.04 | 0.0001 | 7
20 || nan | 0.000 | nan | nan | nan || nan | 0.000 | nan | nan | nan || 1879.82 | 0.0001 | 7
21 || nan | 0.000 | nan | nan | nan || nan | 0.000 | nan | nan | nan || 1878.32 | 0.0001 | 7
22 || nan | 0.000 | nan | nan | nan || nan | 0.000 | nan | nan | nan || 1878.99 | 0.0001 | 7
23 || nan | 0.000 | nan | nan | nan || nan | 0.000 | nan | nan | nan || 1877.28 | 0.0001 | 7
24 || nan | 0.000 | nan | nan | nan || nan | 0.000 | nan | nan | nan || 1881.16 | 0.0000 | 7
25 || nan | 0.000 | nan | nan | nan || nan | 0.000 | nan | nan | nan || 1884.04 | 0.0000 | 7
26 || nan | 0.000 | nan | nan | nan || nan | 0.000 | nan | nan | nan || 1877.73 | 0.0000 | 7
27 || nan | 0.000 | nan | nan | nan || nan | 0.000 | nan | nan | nan || 1882.80 | 0.0000 | 7
No longer model is learning: training stop

处理后可以正常训练，训练输出的 checkpoint 保存在 ./save_models 文件夹中。损失最小的 checkpoint 是，save_32.pt。将复制到新建的 ./best_model，并重命名为 best_model.pt。

3.4 分子生成案例测试

作者提供了两个 jupyter notebook，分别是 DeepICL_Ligand_Elaboration_DEMO.ipynb 和 DeepICL_No_Reference_Ligand_DEMO.ipynb 分别是在有、无参考分子的时候进行分子生成

3.4.1 内置的小分子优化案例

项目提供的小分子优化案例和基于口袋的从头分子生成案例使用的都是 2F0Z 蛋白，原始配体在口袋中的结构如下：

配体的 2D 结构如下：

项目提供基于参考配体的分子优化案例，项目提供示例：DeepICL_Ligand_Elaboration_DEMO.ipynb，在谷歌网盘中，链接为：https://drive.google.com/file/d/1Sg2mIeFut66KjAhZBgM-1nPqSmc9g7lC/view?usp=sharing。创建 ./notebook 文件夹，示例 jupyter notebook 下载后保存到 ./notebook 文件夹中。

当前，整个项目的目录情况如下：

.
|-- LICENSE
|-- README.md
|-- best_model
|-- check_processed_data.ipynb
|-- crossdocked
|-- data
|-- example
|-- image
|-- install_packages.sh
|-- notebook
|-- process.log
|-- processed_data
|-- save_model
|-- script
|-- temp
`-- v2020

11 directories, 5 files

其中，./best_model/save_best.pt 是上述训练好的模型，./notebook 中是下载的 jupyter notebook 示例。./example_elaboration 中是 2f0z 的配体和蛋白结构，分别是 2f0z_ligand.sdf 和 2f0z_protein.pdb，具体如下：

.
|-- 2f0z_ligand.sdf
`-- 2f0z_protein.pdb

2 files

接下来，通过项目提供的 DeepICL_Ligand_Elaboration_DEMO.ipynb 介绍分子优化过程。

导入分子优化所需的依赖库，设置随机种子。

#load dependencies
import os
import numpy as np
import py3Dmol
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
import torch
import torch_geometric
import Bio
import rdkit
import plip
import glob
import random
import multiprocessing

def seed_torch(seed=0):
    random.seed(seed)
    os.environ['PYTHONHASHSEED'] = str(seed)
    np.random.seed(seed)
    torch.manual_seed(seed)
    torch.cuda.manual_seed(seed)
    torch.backends.cudnn.benchmark = False
    torch.backends.cudnn.deterministic = True

seed_torch(2024)

指定参考配体和蛋白质结构的位置，设置生成分子的保存路径和工作路径 dir_name 为 '/workspace/xxxx/projects/DeepICL/example' ，即项目的主目录 ./。

dir_name = '/workspace/xxxx/projects/DeepICL' #
Save_Directory = 'example_elaboration' #@param {type:"string"}
Protein_PDB_file_name = '2f0z_protein.pdb' #@param {type:"string"}
Reference_Ligand_SDF_file_name = '2f0z_ligand.sdf'  #@param {type:"string"}


workDir = f"{dir_name}/{Save_Directory}"
protein_pdb = os.path.join(workDir, str(Protein_PDB_file_name))
ref_ligand_sdf = os.path.join(workDir, str(Reference_Ligand_SDF_file_name))


assert os.path.exists(protein_pdb)
assert os.path.exists(ref_ligand_sdf)

通过 openbabel 把配体和蛋白保存成复合体结构，保存的复合体结构为 /workspace/xxxx/projects/DeepICL/example/tmp_complex.pdb。

!obabel "{ref_ligand_sdf}" "{protein_pdb}" -O ../example_elaboration/tmp_complex.pdb -j

设置使用的 CPU 核心数（NCPU）为 2 、原子类型选择和位置选择的随机因子（temperature_factor_1、temperature_factor_2）。设置使用条件（Use_condition）、分子采样数目设置为 50（Num_sample）等。指定分子骨架为 "C1CCCCO1"。

Result_file_name = "2f0z_results" #@param {type: "string"}
NCPU = "2" #@param {type:"string"}
Temperature_factor = "0.05" #@param {type:"string"}
Use_scaffold = True #@param {type:"boolean"}
Use_condition = True #@param {type:"boolean"}
Num_sample = 50 #@param {type:"slider", min:50, max:500, step:50}
Verbose = True #@param {type:"boolean"}


ncpu = NCPU
result_file_name = f"{workDir}/{Result_file_name}" #
temperature_factor_1 = float(Temperature_factor)
temperature_factor_2 = float(Temperature_factor)
use_scaffold = Use_scaffold
use_condition = Use_condition
num_sample = Num_sample
verbose = Verbose
predefined_scaffold = None
if use_scaffold:
  #@markdown *Optional*
  Scaffold_SMILES = "C1CCCCO1" #@param {type:"string"}
  predefined_scaffold = Scaffold_SMILES

对原始的配体结构和蛋白结构进行数据预处理（参考 3.2 数据预处理），并显示参考配体中指定的分子骨架结构。

#@title **Pre-processing data**


os.chdir(f'{dir_name}/data')
from processor import PDBbindDataProcessor
from rdkit import Chem
import pickle


prep = PDBbindDataProcessor(
    data_dir="tmp",
    save_dir="tmp",
    use_whole_protein=True,
    predefined_scaffold=predefined_scaffold
)
prep._tmp_dir = workDir


data = prep._processor(
    ligand_fn=ref_ligand_sdf,
    pocket_fn=protein_pdb
)


ref_lig_smiles = data[9]
m = Chem.MolFromSmiles(ref_lig_smiles)
if use_scaffold:
    ref_core_smiles = data[10]
    c = Chem.MolFromSmarts(ref_core_smiles)
    indices = m.GetSubstructMatches(c)[0]
else:
    indices = []


Data_save_file_name = "2f0z" #@param {type: "string"}
with open(f"{workDir}/{Data_save_file_name}", 'wb') as w:
    pickle.dump(data, w)
with open(f"{workDir}/test_keys.pkl", 'wb') as w:
    pickle.dump([Data_save_file_name], w)


from rdkit.Chem import Draw
from rdkit.Chem.Draw import rdMolDraw2D
from rdkit.Chem.Draw import IPythonConsole
from IPython.display import SVG, display


drawer = rdMolDraw2D.MolDraw2DSVG(500, 500)
drawer.DrawMolecule(m, highlightAtoms=indices)
drawer.FinishDrawing()
svg = drawer.GetDrawingText()
display(SVG(svg))

运行输出：

3 4 SINGLE

3 12 SINGLE

4 5 SINGLE

5 6 SINGLE

6 11 SINGLE

11 12 SINGLE

基于原始蛋白和配体结构的数据预处理输出文件保存为 ./example_elaboration/2f0z 。保存的 ./example_elaboration/test_keys.pkl 为体系的 PDB ID 索引。

在参考分子中，指定的分子骨架使用红色突出（对应上文中的骨架 C1CCCCO1， 为六元含氧杂环），如下：

设置完整的分子生成命令。

os.chdir(f'../')
command = f"generate.py --ncpu {NCPU} --k 8 --data_dir '{workDir}' --key_dir '{workDir}' " + \
          f"--restart_dir /workspace/xxxx/projects/DeepICL/best_model/save_best.pt --result_dir '{result_file_name}' --num_layers 6 --num_dense_layers 3 --num_hidden_feature 128 --num_cond_feature 7 " + \
          f"--num_sample {num_sample} --max_num_add_atom 30 --dist_one_hot_param1 0 10 25 --dist_one_hot_param2 0 15 300 --temperature_factor1 {temperature_factor_1} " + \
          f"--temperature_factor2 {temperature_factor_2} --radial_limits 0.9 2.2 --pocket_coeff_max 10.0 --pocket_coeff_thr 2.5 --pocket_coeff_beta 0.91 --conditional -y"
if use_scaffold:
    command += " --use_scaffold"
if use_condition:
    command += " --use_condition"
if verbose:
    command += " --verbose"


print(command)

至此已经设置完成分子生成的参数，具体的分子生成命令如下：

python generate.py \
  --ncpu 2 \
  --k 8 \
  --data_dir /workspace/xxxx/projects/DeepICL/example \
  --key_dir /workspace/xxxx/projects/DeepICL/example \
  --restart_dir /workspace/xxxx/projects/DeepICL/best_model/save_best.pt \
  --result_dir /workspace/xxxx/projects/DeepICL/example/2f0z_results \
  --num_layers 6 \
  --num_dense_layers 3 \
  --num_hidden_feature 128 \
  --num_cond_feature 7 \
  --num_sample 1 \
  --max_num_add_atom 30 \
  --dist_one_hot_param1 0 10 25 \
  --dist_one_hot_param2 0 15 300 \
  --temperature_factor1 0.05 \
  --temperature_factor2 0.05 \
  --radial_limits 0.9 2.2 \
  --pocket_coeff_max 10.0 \
  --pocket_coeff_thr 2.5 \
  --pocket_coeff_beta 0.91 \
  --conditional \
  --use_scaffold \
  --use_condition \
  --verbose \
  -y 

--ncpu 2 指定分子生成使用两个 CPU 核心；--k 8 表示使用 8-近邻；--data_dir /workspace/xxxx/projects/DeepICL/example 指定预处理好的数据所在的文件夹；--key_dir /workspace/xxxx/projects/DeepICL/example 指定 test_keys.pkl 所在的文件夹；--restart_dir /workspace/xxxx/projects/DeepICL/best_model/save_best.pt 指定训练好的模型；--result_dir /workspace/xxxx/projects/DeepICL/example/2f0z_results 指定输出结果保存的文件夹；--num_sample 50 表示采样 50 个分子；--temperature_factor1 0.05 和 --temperature_factor2 0.05 分别设置类型选择和位置选择的控制随机性；--conditional 表示在采样时模型使用相互作用作为条件；--use_scaffold 表示使用指定的分子骨架。

运行分子生成命令

os.chdir(f'{dir_name}/script')
print(os.getcwd())
!python {command}

基于 2f0z 的口袋生成 50 个分子，保存在 ./example/2f0z_results 中，花费时间约 2.5 分钟。

/workspace/xxxx/projects/DeepICL/example 文件夹中保存中间处理文件和生成分子，文件结构如下：

.
|-- 2f0z
|-- 2f0z_ligand.sdf
|-- 2f0z_protein.pdb
|-- 2f0z_results
|-- test_keys.pkl
`-- tmp_complex.pdb


1 directories, 5 files

其中，2f0z 是根据原始结构处理的模型的输入信息。2f0z_ligand.sdf 和 2f0z_protein.pdb 分别是配体和蛋白结构。2f0z_results 文件夹中是生成分子，其中每个生成分子都保存为单个的 SDF 文件，还包括一个所有生成分子合并后的 SDF 文件，即 2f0z_novel_results.sdf 。test_keys.pkl 是处理好的 key 。tmp_complex.pdb 是配体和蛋白的复合物结构。

接着，评价生成分子。'result_file_name' 是生成分子的路径，即 ./example_elaboration/2f0z_results。

--key_dir '{workDir}/test_keys.pkl' 指定的测试 pkl 文件位置，为 ./example_elaboration/test_keys.pkl。--smi_dir './data/keys/train_smiles.txt' 指定训练DeepICL的训练集的小分子的smiles（用于计算新颖性）。

!python evaluate.py \
  --result_dir '{result_file_name}' \
  --key_dir '{workDir}/test_keys.pkl' \
  --smi_dir '{dir_name}/data/keys/train_smiles.txt' \
  --filter_level 0

评价生成分子的命令输出：

Key: 2f0z
Generated molecule: 50
Validity: 100.00%
Uniqueness: 100.00%
Novelty: 100.00%
TOTAL VALIDITY: 100.00
TOTAL UNIQUNESS: 100.00
TOTAL NOVELTY: 100.00

评估结果显示，生成的 50 个分子中， 全部生成分子都通过了有效性、唯一性和新颖性过滤，表明生成分子全部是化学有效的分子，并且没有重复生成分子，和训练集分子也没有交集。

计算所有生成分子、参考分子分别和蛋白口袋的相互作用，并计算每个生成分子和参考分子的相互作用相似性，找出相似度最高的生成分子。（相互作用相似性的计算方法，计算每一个生成的小分子与蛋白质的相互作用，然后标记口袋原子每一个原子的相互作用类型）

#@title **Find a ligand with a highest interaction similarity**


from sklearn.metrics.pairwise import cosine_similarity
from tqdm import tqdm


def calc_intr_mat(ligand_fn, protein_fn):
    global prep
    ligand_mol = Chem.SDMolSupplier(ligand_fn)[0]
    protein_mol = Chem.MolFromPDBFile(protein_fn)
    ligand_n = ligand_mol.GetNumAtoms()
    protein_n = protein_mol.GetNumAtoms()
    complex_mol, complex_fn = prep._join_complex(ligand_fn, protein_fn)
    intr_info = prep._get_complex_interaction_info(complex_fn)
    intr_mat = prep._get_pocket_interaction_matrix(ligand_n, protein_n, intr_info)
    intr_mat = intr_mat[:, :-1]
    os.unlink(complex_fn)
    return intr_mat.reshape(1, -1)


def calc_similarity(x1, x2, metric=cosine_similarity):
    return metric(x1, x2)[0][0]


gen_unique_ligand_mols = Chem.SDMolSupplier(
    list(glob.glob(f"{result_file_name}/*_novel_results.sdf"))[0]
)
gen_unique_ligand_fns = [f'{result_file_name}/{m.GetProp("_Name")}.sdf' for m in gen_unique_ligand_mols]
ref_intr_mat = calc_intr_mat(ref_ligand_sdf, protein_pdb)
sim_list = []
for gulf in tqdm(gen_unique_ligand_fns):
    gen_intr_mat = calc_intr_mat(gulf, protein_pdb)
    sim = calc_similarity(ref_intr_mat, gen_intr_mat)
    sim_list.append((sim, gulf))
sim_list = sorted(sim_list, key=lambda x: x[0], reverse=True)


print(f"\n\"{sim_list[0][1].split('/')[-1]}\" showed the highest interaction similarity of {sim_list[0][0]:.3f}.")
gen_ligand_sdf = sim_list[0][1]

计算结果输出：

"2f0z_33.sdf" showed the highest interaction similarity of 0.875.

生成分子中的 2f0z_33.sdf 和参考分子的相互作用相似度最高，为 0.875 。该分子最大程度的保留原始的相互作用。

参考分子 2f0z_ligand 与 相互作用最相似分子 2f0z_33.sdf 的 2D 结构见下图。2D结构具有非常高的相似性。

在口袋空间占据方面，如下图。紫色的是生成的分子 2f0z_33.sdf ，蓝色的是 参考分子2f0z_ligand 。二者的口袋空间占据非常相似，但生成的分子 2f0z_33.sdf 支链更长，占据了更多半开放的口袋空间。

参考分子 2f0z_ligand 与口袋之间的相互作用模式如下左图，相互作用相似度最高的 2f0z_33.sdf 与口袋的相互作用，如右下图。在参考分子中，关键相互作用氨基酸有， Glu39, ARG21, Arg304, Met65, Tyr179。这些氨基酸均出现在了生成的 2f0z_33.sdf 有相互作用氨基酸中，但在延申的支链处形成了新的相互作用（Tyr181, Glu111 ）。

注：为后续测试案例不覆盖生成分子，将 ./example 更名为 ./example_elaboration 。

所有生成的分子如下：

作者没有让生成分子对接到蛋白口袋中，查看在口袋中的结合情况，接下来，我们通过 qvina 对接评估生成分子。在 ./example_elaboration 文件夹中创建 ./example_elaboration/1_get_qvina 文件夹，以存放我们计算 qvina 的脚本和处理文件。

(注：这里略过代码部分)

使用 MOE 打开，并根据 qvina_score 进行排序。qvina score 排名前 3 的分子的 2D 结构如下，对应的 qvina score 分数分别为 -8.8， -8.6， -8.5：

qvina score 排名前 3 的分子在口袋中的 Pose 如下：

3.4.2 内置的从头分子生成案例

项目提供仅有蛋白口袋，没有参考分子的从头分子生成案例，项目提供示例：DeepICL_No_Reference_Ligand_DEMO.ipynb，在谷歌网盘中，链接为：https://drive.google.com/file/d/1Sg2mIeFut66KjAhZBgM-1nPqSmc9g7lC/view 。 下载后保存到 ./notebook 文件夹中。创建 ./example_no_ref，文件夹中包含原始配体和蛋白结构，即 2f0z_ligand.sdf 和 2f0z_protein.pdb 。接下来通过项目提供的示例介绍分子生成过程，以下代码均在 ipynb 文件中执行。

导入相关的依赖库并设置随机种子。

#load dependencies
import os
import numpy as np
import py3Dmol
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
import torch
import torch_geometric
import Bio
import rdkit
import plip
import glob
import random
import multiprocessing


def seed_torch(seed=0):
    random.seed(seed)
    os.environ['PYTHONHASHSEED'] = str(seed)
    np.random.seed(seed)
    torch.manual_seed(seed)
    torch.cuda.manual_seed(seed)
    torch.backends.cudnn.benchmark = False
    torch.backends.cudnn.deterministic = True


seed_torch(2024)

指定输入文件和输出保存路径。

#@title **Please, provide the necessary input files below**:
dir_name = '/workspace/xxx/projects/DeepICL'
Save_Directory = '/example_no_ref' #@param {type:"string"}
Protein_PDB_file_name = '2f0z_protein.pdb' #@param {type:"string"}


workDir = f"{dir_name}/{Save_Directory}"
protein_pdb = os.path.join(workDir, str(Protein_PDB_file_name))


assert os.path.exists(protein_pdb)

设置使用的 CPU 核心数（NCPU）为 2 、指定蛋白质口袋中心坐标（pocket_center_x , pocket_center_y , pocket_center_z）、原子类型选择和位置选择的随机因子（Temperature_factor）。设置使用条件（Use_condition）、分子采样数目设置为 50（Num_sample）、给初始位置添加高斯噪音（Add_noise）等。

Result_file_name = "2f0z_results" #@param {type: "string"}
NCPU = "2" #@param {type:"string"}
Pocket_center_X = 33.86 #@param {type: "string"}
Pocket_center_Y = 26.98 #@param {type: "string"}
Pocket_center_Z = 29.13 #@param {type: "string"}
Temperature_factor = "0.05" #@param {type:"string"}
Use_condition = True #@param {type:"boolean"}
Num_sample = 50 #@param {type:"slider", min:50, max:500, step:50}
Add_noise = True #@param {type:"boolean"}
Verbose = True #@param {type:"boolean"}


ncpu = NCPU
result_file_name = f"{workDir}/{Result_file_name}"
temperature_factor_1 = float(Temperature_factor)
temperature_factor_2 = float(Temperature_factor)
pocket_center_x = float(Pocket_center_X)
pocket_center_y = float(Pocket_center_Y)
pocket_center_z = float(Pocket_center_Z)
pocket_center = [
    pocket_center_x,
    pocket_center_y,
    pocket_center_z
]
use_condition = Use_condition
num_sample = Num_sample
add_noise = Add_noise
verbose = Verbose

处理蛋白口袋结构

os.chdir(f'{dir_name}/data')
from processor_no_ligand import PocketDataProcessor
from rdkit import Chem
import pickle


prep = PocketDataProcessor(
    data_dir="tmp",
    save_dir="tmp",
)
prep._tmp_dir = workDir


data = prep._processor(
    pocket_fn=protein_pdb,
    pocket_center=pocket_center
)


Data_save_file_name = "2f0z_no_ligand" #@param {type: "string"}
with open(f"{workDir}/{Data_save_file_name}", 'wb') as w:
    pickle.dump(data, w)
with open(f"{workDir}/test_keys.pkl", 'wb') as w:
    pickle.dump([Data_save_file_name], w)

设置完整的分子生成命令

os.chdir(f'../')
command = f"./generate.py --ncpu {NCPU} --k 8 --data_dir '{workDir}' --key_dir '{workDir}' " + \
          f"--restart_dir '{dir_name}'/best_model/save_best.pt --result_dir '{result_file_name}' --num_layers 6 --num_dense_layers 3 --num_hidden_feature 128 --num_cond_feature 7 " + \
          f"--num_sample {num_sample} --max_num_add_atom 30 --dist_one_hot_param1 0 10 25 --dist_one_hot_param2 0 15 300 --temperature_factor1 {temperature_factor_1} " + \
          f"--temperature_factor2 {temperature_factor_2} --radial_limits 0.9 2.2 --pocket_coeff_max 10.0 --pocket_coeff_thr 2.5 --pocket_coeff_beta 0.91 --conditional -y"
if use_condition:
    command += " --use_condition"
if add_noise:
    command += " --add_noise"
if verbose:
    command += " --verbose"


print(command)

至此已经设置完成分子生成的参数，具体的分子生成命令如下：

python ./generate.py \
  --ncpu 2 \
  --k 8 \
  --data_dir /workspace/xxxx/projects/DeepICL/example \
  --key_dir /workspace/xxxx/projects/DeepICL/example \
  --restart_dir /workspace/xxxx/projects/DeepICL/best_model/save_best.pt \
  --result_dir /workspace/xxxx/projects/DeepICL/example/2f0z_results \
  --num_layers 6 \
  --num_dense_layers 3 \
  --num_hidden_feature 128 \
  --num_cond_feature 7 \
  --num_sample 50 \
  --max_num_add_atom 30 \
  --dist_one_hot_param1 0 10 25 \
  --dist_one_hot_param2 0 15 300 \
  --temperature_factor1 0.05 \
  --temperature_factor2 0.05 \
  --radial_limits 0.9 2.2 \
  --pocket_coeff_max 10.0 \
  --pocket_coeff_thr 2.5 \
  --pocket_coeff_beta 0.91 \
  --conditional \
  --use_condition \
  --add_noise \
  --verbose \
  -y 

--ncpu 2 指定分子生成使用两个 CPU 核心；--k 8 表示使用 8-近邻；--data_dir /workspace/xxxx/projects/DeepICL/example 指定预处理好的数据所在的文件夹；--key_dir /workspace/xxxx/projects/DeepICL/example 指定 test_keys.pkl 所在的文件夹；--restart_dir /workspace/xxxx/projects/DeepICL/best_model/save_best.pt 指定训练好的模型；--result_dir /workspace/xxxx/projects/DeepICL/example/2f0z_results 指定输出结果保存的文件夹；--num_sample 50 表示采样 50 个分子。

--temperature_factor1 0.05 和 --temperature_factor2 0.05 分别设置类型选择和位置选择的控制随机性；--conditional 表示在采样时模型使用相互作用作为条件；--add_noise 表示给初始的位置添加高斯噪音。

运行分子生成命令

os.chdir(f'{dir_name}/script')
print(os.getcwd())
!python {command}

基于 2f0z 的口袋生成 50 个分子，保存在 ./example/2f0z_results 中，花费时间约 2 分钟。/workspace/xxxx/projects/DeepICL/example 文件夹中保存中间处理文件和生成分子，文件结构如下：

.
|-- 2f0z_ligand.sdf
|-- 2f0z_no_ligand
|-- 2f0z_on_ref_results_all.sdf
|-- 2f0z_protein.pdb
|-- 2f0z_protein.pdbqt
|-- 2f0z_results
`-- test_keys.pkl


1 directories, 6 files

其中，2f0z_ligand.sdf 和 2f0z_protein.pdb 分别是配体和蛋白结构。2f0z_results 文件夹中是生成分子，其中每个生成分子都保存为单个的 SDF 文件，还包括一个所有生成分子合并后的 SDF 文件，即 2f0z_no_ligand_novel_results.sdf 。test_keys.pkl 是处理好的 key 。

接着，评价生成分子。'result_file_name' 是生成分子的路径，即 /workspace/xxxx/projects/DeepICL/example/2f0z_results。

--key_dir '{workDir}/test_keys.pkl' 指定的测试 pkl 文件位置，为 /workspace/xxxx/projects/DeepICL/example/test_keys.pkl。--smi_dir '{dir_name}/data/keys/train_smiles.txt' 指定训练 DeepICL 模型时训练集的小分子smiles，保存在 /workspace/xxxx/projects/DeepICL/data/keys/test_smiles.txt 。

!python evaluate.py \
  --result_dir '{result_file_name}' \
  --key_dir '{workDir}/test_keys.pkl' \
  --smi_dir '{dir_name}/data/keys/train_smiles.txt' \
  --filter_level 0

评价生成分子的命令输出：

Key: 2f0z_no_ligand
Generated molecule: 50
Validity: 98.00%
Uniqueness: 100.00%
Novelty: 100.00%
TOTAL VALIDITY: 98.00
TOTAL UNIQUNESS: 100.00
TOTAL NOVELTY: 100.00

评估结果显示，生成的 50 个分子中，49 个分子通过有效性过滤，分子有效性达 98% 。 全部生成分子都通过了唯一性和新颖性过滤，表明生成分子中没有重复分子，和训练集分子也没有交集。

所有生成的分子如下：

作者没有让生成分子对接到蛋白口袋中，查看在口袋中的结合情况，接下来，我们通过 qvina 对接评估生成分子。在 ./example_no_ref 文件夹中创建 ./example_no_ref/1_get_qvina 文件夹，以存放我们计算 qvina 的脚本和处理文件。（注：略过，代码相见附件资源）

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

qvina score 排名前 3 的分子在口袋中的 Pose 如下：

我们在 ./example_no_ref 文件夹中创建 ./example_no_ref/get_highest_interaction_similarity.ipynb 来计算每个生成分子和参考分子分别与口袋的相互作用，并挑选出和参考分子相互作用相似度最高的生成分子。具体代码如下：

导入相关的依赖库，设置生成分子的文件路径，即 ./example_no_ref /2f0z_results。配体结构和蛋白结构的位置分别为：/workspace/xxxx/projects/DeepICL/example_no_ref/2f0z_ligand.sdf 和 /workspace/xxxx/projects/DeepICL/example_no_ref/2f0z_protein.pdb 。

from rdkit.Chem import Draw
import os


os.chdir(f'/workspace/xxxx/projects/DeepICL/data')
from processor import PDBbindDataProcessor
from rdkit import Chem
import pickle


import numpy as np
import py3Dmol
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
import torch
import torch_geometric
import Bio
import rdkit
import plip
import glob
import random
import multiprocessing


result_file_name = '/workspace/xxxx/projects/DeepICL/example_no_ref/2f0z_results'


dir_name = '/workspace/xxxx/projects/DeepICL' # xxxx
Save_Directory = 'example_no_ref' #@param {type:"string"}
Protein_PDB_file_name = '2f0z_protein.pdb' #@param {type:"string"}
Reference_Ligand_SDF_file_name = '2f0z_ligand.sdf'  #@param {type:"string"}


workDir = f"{dir_name}/{Save_Directory}"
protein_pdb = os.path.join(workDir, str(Protein_PDB_file_name))
ref_ligand_sdf = os.path.join(workDir, str(Reference_Ligand_SDF_file_name))

计算每个生成分子和蛋白口袋的相互作用向量，原始配体和蛋白口袋的相互作用向量。计算所有生成分子相互作用向量和参考分子的余弦相似度，挑选出相似度最高的生成分子。

prep = PDBbindDataProcessor(
    data_dir="tmp",
    save_dir="tmp",
    use_whole_protein=True
)
prep._tmp_dir = workDir


#@title **Find a ligand with a highest interaction similarity**


from sklearn.metrics.pairwise import cosine_similarity
from tqdm import tqdm


def calc_intr_mat(ligand_fn, protein_fn):
    global prep
    ligand_mol = Chem.SDMolSupplier(ligand_fn)[0]
    protein_mol = Chem.MolFromPDBFile(protein_fn)
    ligand_n = ligand_mol.GetNumAtoms()
    protein_n = protein_mol.GetNumAtoms()
    complex_mol, complex_fn = prep._join_complex(ligand_fn, protein_fn) # 返回的是复合物的 mol 对象和 文件名
    intr_info = prep._get_complex_interaction_info(complex_fn) # (anion_indices,cation_indices,hbd_indices,hba_indices,hyd_indices,pipi_indices,None,) 返回的是一个元组，包含了复合物的相互作用信息
    intr_mat = prep._get_pocket_interaction_matrix(ligand_n, protein_n, intr_info) # 返回的是一个矩阵，矩阵的每一行代表了一个原子和蛋白质的相互作用信息
    intr_mat = intr_mat[:, :-1] # 把最后一列去掉，因为最后一列是 None
    os.unlink(complex_fn) # 删除复合物文件
    return intr_mat.reshape(1, -1) # 把矩阵展平


def calc_similarity(x1, x2, metric=cosine_similarity):
    return metric(x1, x2)[0][0]


gen_unique_ligand_mols = Chem.SDMolSupplier(
    list(glob.glob(f"{result_file_name}/*_novel_results.sdf"))[0]
)
gen_unique_ligand_fns = [f'{result_file_name}/{m.GetProp("_Name")}.sdf' for m in gen_unique_ligand_mols] # 生成的小分子的文件名列表
ref_intr_mat = calc_intr_mat(ref_ligand_sdf, protein_pdb) # 参考小分子和蛋白质的相互作用矩阵，这个矩阵是一个行向量


# 遍历生成分子的文件名列表，计算每个生成分子和参考分子的相互作用余弦相似度
sim_list = []
for gulf in tqdm(gen_unique_ligand_fns):
    gen_intr_mat = calc_intr_mat(gulf, protein_pdb)
    sim = calc_similarity(ref_intr_mat, gen_intr_mat)
    sim_list.append((sim, gulf))
sim_list = sorted(sim_list, key=lambda x: x[0], reverse=True)


print(f"\n\"{sim_list[0][1].split('/')[-1]}\" showed the highest interaction similarity of {sim_list[0][0]:.3f}.")
gen_ligand_sdf = sim_list[0][1]

输出结果：

"2f0z_no_ligand_2.sdf" showed the highest interaction similarity of 0.645.

生成分子中的 2f0z_no_ligand_2.sdf 和参考分子的相互作用相似度最高，为 0.645 。该分子最大程度的保留原始的相互作用。参考分子 2f0z_ligand 与 相互作用最相似分子 2f0z_no_ligand_2.sdf 的 2D 结构见下图。

在口袋空间占据方面，如下图。蓝色的是生成的分子 2f0z_no_ligand_2.sdf ，绿色的是 参考分子2f0z_ligand 。二者的口袋空间中心都是一个环结构占据，但生成的分子 2f0z_no_ligand_2.sdf 具有更多朝向口袋外的支链，占据了更多半开放的口袋空间。

参考分子 2f0z_ligand 与口袋之间的相互作用模式如下左图，相互作用相似度最高的 2f0z_no_ligand_2.sdf 与口袋的相互作用，如右下图。在参考分子中，关键相互作用氨基酸有， Glu39, ARG21, Arg304, Met85, Tyr179。但 2f0z_no_ligand_2.sdf 和口袋形成相互作用的氨基酸只有 Glu111、Glu218 和 Arg304 。从头分子生成并没有参考分子，只会根据口袋的形状和环境来生成分子，所以和参考分子的相互作用相似性没有办法保证。但是会生成更多新颖的分子。

3.4.3 自定义的从头分子生成案例（no_reference）

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

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

具体的生成过程，这里略过，相见附件资源。这里直接讲结果。

DeepICL生成了50个分子， 评价生成分子的命令输出：

Key: 3wze_no_ligand
Generated molecule: 50
Validity: 96.00%
Uniqueness: 100.00%
Novelty: 100.00%
TOTAL VALIDITY: 96.00
TOTAL UNIQUNESS: 100.00
TOTAL NOVELTY: 100.00

评估结果显示，生成的 50 个分子中，48 个分子通过有效性过滤，分子有效性达 96% 。 全部生成分子都通过了唯一性和新颖性过滤，表明生成分子中没有重复分子，和训练集分子也没有交集。

所有生成的分子如下：

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

qvina score 排名前 3 的分子在口袋中的 Pose 如下：

生成分子中的 3wze_no_ligand_25.sdf 和参考分子的相互作用相似度最高，为 0.669 。该分子最大程度的保留原始的相互作用。参考分子 3wze_ligand 与 相互作用最相似分子 3wze_no_ligand_25.sdf 的 2D 结构见下图。

在口袋空间占据方面，如下图。蓝色的是生成的分子 3wze_no_ligand_25.sdf ，绿色的是 参考分子 3wze_ligand 。原始的参考配体的链相较于 3wze_no_ligand_25.sdf 更长，能够占据右端口袋。

参考分子 3wze_ligand 与口袋之间的相互作用模式如下左图，相互作用相似度最高的 3wze_no_ligand_25.sdf 与口袋的相互作用，如右下图。在参考分子中，关键相互作用氨基酸有：Cys919、Leu1035、Glu885、IIe1044、Cys1045、Asp1046、Phe1047 。3wze_no_ligand_25.sdf 和口袋形成相互作用的氨基酸：Asp1046、Cys919、Phe918 。我们自定义的测试案例，没有指定参考分子，只根据蛋白口袋生成分子。所以在保留参考配体的相互作用方面表现较差。

3.4.4 自定义的从头分子生成案例（reference）

使用类似 3.4.2 的方法，可以生成 3wze 体系的小分子。首先定义分子的参考骨架是（下图红色区域）：

评估结果显示，生成的 50 个分子中，50 个分子通过有效性过滤，分子有效性达 100% 。 全部生成分子都通过了唯一性和新颖性过滤，表明生成分子中没有重复分子，和训练集分子也没有交集。

Key: 3wze
# Generated molecule: 50
Validity: 100.00%
Uniqueness: 100.00%
Novelty: 100.00%
TOTAL VALIDITY: 100.00
TOTAL UNIQUNESS: 100.00
TOTAL NOVELTY: 100.00

相互作用模式相似性检测，最高相似为0.933，对应 3wze_16.sdf 小分子。

"3wze_16.sdf" showed the highest interaction similarity of 0.933

如下图，绿色为生成的分子，紫色为参考分子，二者药效团分布还是非常的像，形状之间相似度也很高：

3wze_16.sdf 与蛋白的相互作用如下，除了core 部分，与参考分子的相互作用没有非常相似，在原来的 Cys191 附近仍保持了一个 与 Gly922 的氢键。

所有生成分子经过 qvina 对接以后，最高打分的三个分子的打分分别为：-13.3, -12.2, -11.8。对应的分子在口袋中的形状如下：

所有生成的分子如下：

从这些生成分子的情况来看，生成的小分子对口袋空间的占据比较充分，同时药效团原子与参考分子的比较相似，相互作用模式大部分得以保留，生成分子的 qvina 打分也不错。但是生成分子的构象存在一定的问题，不饱和环，特殊环环状结构的存在。 

四、总结

DeepICL 是一种相互作用感知的自回归的基于等变神经网络的 3D 分子生成 条件 VAE 框架，该框架整合了蛋白质-配体相互作用的先验知识，以泛化基于结构的药物设计。与之前仅依赖于结合结构信息的模型相比，作者更注重表征四种类型的蛋白质-配体相互作用的泛化模式（疏水相互作用、氢键、盐桥和 π-π 堆积）。作者通过全面分析为未见过的目标设计的配体在各个方面的表现，展示了方法的泛化能力，并确认框架能够以高比例建立有利的相互作用，这得益于条件生成框架的可控制性。

根据我们上述的测试结果，在有参考分子的情况下，生成的小分子与参考分子有较多相似的地方，相互作用模式也有一定相似，一定程度上相互作用条件的分子生成的可行性。但是生成的分子结构及其构象还是非常的不理想。

完整的测评文档及此项目可运行的代码详见附件资源。