做复印机的模板网站,wordpress前台,青岛建站程序,制作网站费用分类DeepICL 是一个基于相互作用感知的 3D 分子生成模型#xff0c;能够在目标结合口袋内进行相互作用引导的小分子设计。DeepICL 通过利用蛋白质-配体相互作用的普遍模式作为先验知识#xff0c;在有限的实验数据下也能实现高度的泛化能力。 一、背景介绍 DeepICL 来源于韩国科学…
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 基于相互作用条件的配体设计效果
作者首先展示了相互作用感知条件在设计特定相互作用模式上的效果。参考蛋白质-配体化合物相互作用模式来评估生成分子在相互作用模式方面的表现。作者选择了三种蛋白质分别是骨形态发生蛋白1BMP1、成纤维细胞生长因子1FGF1和二氢叶酸还原酶DHFR。从它们原始复合物的参考配体中提取了初始核心结构。核心结构是参考配体去除了链和官能团保留了由单环或双环组成的最小结构如下图 a 所示。从下图 b 可以看出设计的新配体和参考配体具有很高的相互作用相似性。 下图 c 和 d 分别展示了原始配体和设计的配体具有最高相互作用相似性的相互作用图。对于第一个案例模型成功设计了氢键、π–π 堆积和盐桥符合给定的相互作用条件但具有不同的结构特征。值得注意的是下图 d-1 显示模型能够生成噻吩环代替原始的苯环以与酪氨酸 68TYR68形成 π–π 堆积。这表明模型学习到了芳香族基团所需的特征以形成 π–π 堆积。虽然模型在疏水的苯丙氨酸 157PHE157附近添加了脂肪族碳但距离略大于被概况为疏水相互作用的阈值。对于其他两个案例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 python3.7.13 -y
conda activate DeepICL# 2. install packages dependent on mkl
conda install scipy1.6.2 numpy1.21.2 pandas1.3.4 scikit-learn1.0.2 seaborn0.11.0 -y# 3. install pytorch and torch_geometric
conda install pytorch1.9.0 cudatoolkit10.2 -c pytorch -y
conda install pyg2.0.3 -c pyg -c conda-forge -y# 4. install others
conda install -c rdkit rdkit2022.09.1 -y
conda install -c conda-forge biopython1.77 openbabel3.1.1 -y
conda install -c conda-forge plip2.2.2 -y3.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
-- v20209 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 _processorassert 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 _processorassert self._filter_ligand(ligand_mol), ligand_fn.split(/)[-1]
AssertionError: 6pek_ligand.sdf6pek
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_nsmilesligand_smi相互作用信息作为口袋的条件pocket_cond配体质心center_of_mass 。下面查看主要的内容。 data[0] 为配体的原子类型ligand_typeATOM_TYPES [C, N, O, F, P, S, Cl, Br] 和一个未知项9 维的 one-hot 向量表示一个原子的类型ligand 包含 23 个非氢原子矩阵形状为 239。具体内容如下
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矩阵形状为 233。内容如下
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 个非氢原子矩阵形状为 16251。数据示意如下 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矩阵形状为 1623。数据示意如下
array([[-4.17738318e-01, -8.11868224e00, -2.58775701e00],[-1.42973832e00, -7.51268224e00, -3.42775701e00],[-1.20773832e00, -7.82368224e00, -4.90775701e00],...[ 2.26261682e-01, 2.00331776e00, -8.69775701e00],[ 2.15826168e00, 3.24331776e00, -9.47275701e00],[ 1.59926168e00, 2.08031776e00, -8.93875701e00]]) data[4] 对应配体的邻接矩阵ligand_adj矩阵形状为 2323。数据示意如下
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]], dtypeint32)data[7] 和 data[9] 分别对应配体的原子数目和 SMILES即 23 和 CCCC[CH](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 向量表示。相互作用矩阵的形状为1627;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_parametersTrue 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?uspsharing。创建 ./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
-- v202011 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.pdb2 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 multiprocessingdef seed_torch(seed0):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 Falsetorch.backends.cudnn.deterministic Trueseed_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、分子采样数目设置为 50Num_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 pickleprep PDBbindDataProcessor(data_dirtmp,save_dirtmp,use_whole_proteinTrue,predefined_scaffoldpredefined_scaffold
)
prep._tmp_dir workDirdata prep._processor(ligand_fnref_ligand_sdf,pocket_fnprotein_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, displaydrawer rdMolDraw2D.MolDraw2DSVG(500, 500)
drawer.DrawMolecule(m, highlightAtomsindices)
drawer.FinishDrawing()
svg drawer.GetDrawingText()
display(SVG(svg))
运行输出
3 4 SINGLE3 12 SINGLE4 5 SINGLE5 6 SINGLE6 11 SINGLE11 12 SINGLE 基于原始蛋白和配体结构的数据预处理输出文件保存为 ./example_elaboration/2f0z 。保存的 ./example_elaboration/test_keys.pkl 为体系的 PDB ID 索引。
在参考分子中指定的分子骨架使用红色突出对应上文中的骨架 C1CCCCO1 为六元含氧杂环如下 设置完整的分子生成命令。
os.chdir(f../)
command fgenerate.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 --verboseprint(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.pdb1 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 tqdmdef calc_intr_mat(ligand_fn, protein_fn):global prepligand_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, metriccosine_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, keylambda x: x[0], reverseTrue)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 。 所有生成的分子如下 DeepICL_2f0z_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 multiprocessingdef seed_torch(seed0):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 Falsetorch.backends.cudnn.deterministic Trueseed_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、分子采样数目设置为 50Num_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 pickleprep PocketDataProcessor(data_dirtmp,save_dirtmp,
)
prep._tmp_dir workDirdata prep._processor(pocket_fnprotein_pdb,pocket_centerpocket_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 --verboseprint(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.pkl1 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% 。 全部生成分子都通过了唯一性和新颖性过滤表明生成分子中没有重复分子和训练集分子也没有交集。 所有生成的分子如下 DeepICL_2f0z_no_ref 作者没有让生成分子对接到蛋白口袋中查看在口袋中的结合情况接下来我们通过 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 osos.chdir(f/workspace/xxxx/projects/DeepICL/data)
from processor import PDBbindDataProcessor
from rdkit import Chem
import pickleimport 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 multiprocessingresult_file_name /workspace/xxxx/projects/DeepICL/example_no_ref/2f0z_resultsdir_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_dirtmp,save_dirtmp,use_whole_proteinTrue
)
prep._tmp_dir workDir#title **Find a ligand with a highest interaction similarity**from sklearn.metrics.pairwise import cosine_similarity
from tqdm import tqdmdef calc_intr_mat(ligand_fn, protein_fn):global prepligand_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] # 把最后一列去掉因为最后一列是 Noneos.unlink(complex_fn) # 删除复合物文件return intr_mat.reshape(1, -1) # 把矩阵展平def calc_similarity(x1, x2, metriccosine_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, keylambda x: x[0], reverseTrue)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% 。 全部生成分子都通过了唯一性和新颖性过滤表明生成分子中没有重复分子和训练集分子也没有交集。 所有生成的分子如下 DeepICL-3wze-no-ref 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。对应的分子在口袋中的形状如下 所有生成的分子如下 DeepICL-3wze-all 从这些生成分子的情况来看生成的小分子对口袋空间的占据比较充分同时药效团原子与参考分子的比较相似相互作用模式大部分得以保留生成分子的 qvina 打分也不错。但是生成分子的构象存在一定的问题不饱和环特殊环环状结构的存在。 四、总结
DeepICL 是一种相互作用感知的自回归的基于等变神经网络的 3D 分子生成 条件 VAE 框架该框架整合了蛋白质-配体相互作用的先验知识以泛化基于结构的药物设计。与之前仅依赖于结合结构信息的模型相比作者更注重表征四种类型的蛋白质-配体相互作用的泛化模式疏水相互作用、氢键、盐桥和 π-π 堆积。作者通过全面分析为未见过的目标设计的配体在各个方面的表现展示了方法的泛化能力并确认框架能够以高比例建立有利的相互作用这得益于条件生成框架的可控制性。 根据我们上述的测试结果在有参考分子的情况下生成的小分子与参考分子有较多相似的地方相互作用模式也有一定相似一定程度上相互作用条件的分子生成的可行性。但是生成的分子结构及其构象还是非常的不理想。 完整的测评文档及此项目可运行的代码详见附件资源。
