相互作用感知的蛋白-小分子对接模型

Interformer 是一个应用于分子对接和亲和力预测的深度学习模型，基于 Graph-Transdormer 架构的模型，利用相互作用（氢键、疏水）感知的混合密度网络（interaction-aware mixture den sity network， MDN），捕捉蛋白和小分子之间的非共价相互作用。此外，Interformer 引入负采样策略，有利于有效修正相互作用分布以进行亲和力预测。在PoseBuster 数据集上预测的 top1 构象与晶体构象 RMSD 小于 2A 的比例有 84.09%，在 PDBbind 时间划分的数据集上达到了 63.9%（略差与 SurfDock，相同的数据集，SurfDock 可以达到68.4%）。根据作者描述，在内部的管线上，曾识别到活性分别为 0.7 和16 nM的小分子。

一、背景介绍

Interformer 主要由腾讯深圳的 AI Lab 团队开发，发表文章为《Interformer: an interaction-aware model for protein-ligand docking and affinity prediction》。文章链接：Interformer: an interaction-aware model for protein-ligand docking and affinity prediction | Nature Communications 。该文章在 2024 年 11 月 25 日发表于《Nature Communications》上。

许多蛋白质-配体对接和亲和力预测的深度学习模型忽视了配体与蛋白质原子之间相互作用的复杂建模，导致它们在泛化能力和可解释性方面存在局限性。

在本研究中，作者提出了 Interformer，一种基于 Graph-Transformer 架构、关注相互作用的蛋白质-配体对接与亲和力预测模型。该模型旨在通过利用交互感知混合密度网络（MDN）来捕捉非共价相互作用。此外，作者还引入了一种负采样策略，有效地修正了亲和力预测中的相互作用分布。作者在广泛使用的公开数据集和自有数据集上的实验结果证明了该方法的有效性和普适性。大量分析验证了作者的观点，即通过准确建模特定的蛋白质-配体相互作用，作者的方法能够提升性能。作者的方法在对接任务中达到了最新的（SOTA）性能。

二、模型介绍

近年来，深度学习（DL）方法在分子建模中的应用引起了广泛关注，将对接视为生成建模问题。比如 DiffDock，这是一个基于图神经网络（GNN）的模型，并在结合构象生成方面设立了基准。然而，现有的深度学习模型往往忽视了蛋白质与配体原子之间非共价相互作用的建模，而这些相互作用对于模型的可解释性和泛化能力至关重要。如下图中左图所示，DiffDock 生成的对接构象虽然与晶体结构高度相似，但未能捕捉到非共价相互作用。此外，尽管传统的亲和力预测方法在晶体结构中表现出色，但当面对不太精确的结合构象时，它们的性能急剧下降，这对实际应用提出了挑战。

面对上述的挑战，作者提出了 Interformer。这是一个 AI 模型，旨在解决蛋白质-配体对接中的交互感知问题，并采用构建性学习进行亲和力预测，以应对实际应用中的挑战。

首先，作者提出了一种交互感知混合密度网络（MDN），用于建模非共价相互作用，重点关注蛋白质-配体晶体结构中的氢键和疏水相互作用。

其次，作者引入了伪 Huber 损失函数，利用对比学习指导模型区分有利和不利的结合构象。

第三，所提出的模型基于 Graph-Transformer 框架，已在各种图表示学习任务中证明其优越性，超过了基于 GNN 的其他模型。Interformer 的另一个优势是通过检查 MDN 的融合系数，解释蛋白质-配体相互作用的内部机制。

当在Posebusters 基准数据集和 PDBbind 时间划分基准数据集上评估时，Interformer 在 Posebusters 基准数据集上达到了 84.09% 的准确率，在 PDBbind 时间划分基准数据集上取得了 63.9% 的准确率（以 RMSD 小于 2Å 作为准确的标准）。这一改进归功于模型在捕捉配体和蛋白质之间非共价相互作用方面的增强能力，这对于生成不那么模糊的构象至关重要，并对下游任务的成功表现起着决定性作用。此外，Interformer 甚至在结合构象不够精确的情况下，也能预测合理的亲和力值。作者在自有的实际基准数据集上的评估结果表明，Interformer 与其他模型表现相当，验证了其对构象敏感且具有稳健泛化能力。通过将其应用于实际的内部药物管线中，作者成功地识别了两个小分子，在各自的项目中，它们的亲和力 IC50 值分别为 0.7 nM 和 16 nM，从而证明了其在推动治疗开发中的实际价值。

2.1 模型框架

Interformer 是一个判别式回归模型，基于晶体结构数据进行蛋白质-配体对接任务的训练，并在亲和力预测任务中重新对接对接构象得到对应的亲和力值。核心架构灵感来自 Graph-Transformer，最初用于图表示学习任务。

作者对 Graph-Transformer 进行了修改，引入 Maskedself-attention (MSA) 机制，限制自注意力机制的范围，并在此基础上实现内部（分子内部或蛋白内部）注意力更新机制 Intra-Blocks、外部（蛋白-小分子）更新机制 Inter-Blocks。

在蛋白-小分子的对接任务中，输出边的特征，以预测蛋白-小分子原子对能量。

基于 Graph-Transformer 架构，Interformer 模型管线分为三个阶段：

第一阶段，以单一的初始配体 3D 构象和蛋白质结合位点作为输入，这些数据来自晶体结构。图在各种方法中广泛用于表示配体和蛋白质，如下图所示，其中节点代表原子，边表示两个原子之间的接近关系。作者采用药效团原子类型作为节点特征（22种药效团特征），并使用两个原子之间的欧氏距离作为边特征。这些药效团原子类型提供了重要的化学信息，使模型能够更好地理解特定的相互作用，例如氢键或疏水相互作用。

在第二阶段，如下图所示，处理来自蛋白质和配体的节点特征和边特征，经过 Intra-Blocks 处理。Intra-Blocks 旨在通过捕捉同一分子内的内相互作用来更新每个原子的节点特征。然后，这些更新后的节点特征会输入到 InterBlocks 中，后者旨在捕捉蛋白质和配体原子对之间的相互作用，从而进一步更新节点和边特征。接着，边输出层将这两组特征结合起来，生成每个蛋白质-配体原子对的 Inter 表示。

随后，Inter 表示会通过交互感知的 MDN 进行处理。该网络为每个蛋白质-配体原子对预测四个高斯函数的参数，每个高斯函数分别约束不同类型的特定相互作用。前两个高斯函数封装了所有类型的配对相互作用，第三个高斯函数表示疏水相互作用，第四个高斯函数表示氢键相互作用。通过整合这四个高斯函数，推导出一个混合密度函数（MDF），表示任何给定蛋白质-配体原子对的距离条件概率密度函数。这个 MDF 可以作为能量函数，用来估算蛋白质原子与其对应配体原子之间最可能的距离。

最后，所有蛋白质-配体对的 MDF 会聚合成能量函数的总和，并将其引入到蒙特卡洛（MC）采样方法中，以生成与其蛋白质靶标相关的 top-k 候选配体构象。MC 采样首先将配体放置在蛋白质结合位点的不同位置，并为其分配随机的扭转角度，然后尝试最小化给定的能量函数，以优化配体构象。通过聚合所有根据能量值排序的候选构象，可以获得 top-k 候选对接构象。

在第三阶段，对接得分和亲和力预测管道如下图所示。从生成的对接构象中，蛋白质和配体原子之间的距离和特定相互作用更新新的边特征。然后，节点和边特征会经过 Intra 和 Inter-Blocks 处理，以创建隐性相互作用。一个虚拟节点通过自注意力机制收集所有有关结合构象的信息。最终，虚拟节点的结合嵌入被输入到亲和力和对接得分层中，预测对应对接构象的结合亲和力值和置信度得分。通过引入不良构象，使用对比伪 Huber 损失函数来指导模型判断一个构象是好是坏。训练目标确保模型对不良构象预测较低的值，对良好构象预测较高的值。良好构象与不良构象的主要区别在于它们的相互作用。这一策略帮助模型学习关键的相互作用，而不是人工特征。作者将这一特性称为“对接敏感性”。

2.2 模型性能

下表展示了 Interformer 方法在两个场景下显著超越所有之前的方法，达到了63.9% 的 top-1 成功率，显著高于现有的 SOTA 方法 DiffDock 和 GNINA 。在包括对接得分模型后，top-1 成功率略微下降至 62.1% 。

在PoseBusters 基准数据集中，Interformer 显著优于各种 SOTA AI 和传统模型，成功率达到了 84.09% 。然而，7.8% 的生成构象未通过 PoseBusters 的有效性检查，主要是由于蛋白质和配体原子之间的立体冲突。

在 PDBBind 时间划分测试集中根据序列相似性划分三个子集，分别代表低、中和高同源性水平。作者在这些子集内评估了对接准确性，结果如下图所示。

2.3 特定相互作用复现

在 PDBBind 时间划分测试集的评估中，DiffDock 和 DeepDock 仅能恢复平均 29.42% 、23.55% 的氢键和 19.36% 、16.26% 的疏水相互作用。相比之下，带有对接得分的 Interformer 能恢复平均 57.25% 的氢键和 43.7% 的疏水相互作用。然而，如果没有对接得分，氢键和疏水相互作用的平均恢复率分别略微下降至 52.7% 和 41.6% 。

2.4 负样本提升亲和力预测的姿势敏感性

作者评估了两种不同的亲和力预测策略：一种是包含负样本，另一种仅使用晶体结构而不包含负样本。在 PDBBind 时间划分测试集上，仅使用晶体结构的亲和力模型在预测亲和力值与 RMSD 之间的 Pearson 相关系数为 R=-0.174，说明亲和力预测结果不代表对接 pose 的正确性。然而，当使用负样本时（Affinity w/o Neg），亲和力模型的相关系数为 R=-0.562 ，且带有对接得分模型的相关系数为 R=-0.659 。此外，下图展示了亲和力模型在没有负样本时对好的和差的结合构象做出的一致预测，因为该模型没有利用任何非共价相互作用特征。相反，当引入负样本时，亲和力和对接得分模型都会为具有较大 RMSD 值的结合构象预测较低的亲和力值。结果展示了 Interformer 区分好的和差的构象的能力。

2.5 在真实世界测试集上的亲和力预测评估

作者进一步评估了 Interformer 在 ChEMBL-Kinase， LSD1 project， Mpro covalent test， Mpro project 上的能力，由于姿势敏感性，Interformer 模型在泛化能力上显著优于仅在晶体结构上训练的模型，也优于传统的模型/方法。

为了验证 Interformer 在真实世界场景中的有效性，并展示团队的药物开发能力，作者独立开发了两个药物开发管道。这两个项目都涉及小分子优化，其中药物化学专家基于参考小分子的晶体结构和结合模式设计了一系列候选小分子。随后，Interformer 对这些候选分子进行亲和力评分，并通过各种 ADMET（吸收、分布、代谢、排泄和毒性）分子属性预测模型，确保最终设计的小分子在属性和亲和力方面都表现良好。活性值与预测值之间的关系见下图：

三、Interformer 评测

3.1 安装环境

复制代码项目：

git clone https://github.com/tencent-ailab/Interformer.git

项目提供安装环境的 environment.yml 文件，通过下面命令安装并激活 Interformer 环境：

cd Interformer
conda env create --file environment.yml
conda activate interformer

手动安装其他依赖库，命令如下：

# Install PLIP without wrong openbabel dependency
pip install --no-deps plip# Compile Docking sampling program from source code
cd docking && python setup.py install && cd ../

项目提供训练好的模型和测试数据集等，保存在 zenodo 网站，链接为：Interformer: An Interaction-Aware Model for Protein-Ligand Docking and Affinity Prediction

具体内容截图如下，下载解压后并上传到项目目录（./）。

3.2 对接构象生成案例测试

3.2.1 内置案例

项目提供内置的测试案例是 2qbr.pdb 蛋白，配体在口袋中的构象如下：

配体的 2D 结构如下：

项目提供模型输入的配体和蛋白原始结构文件，保存在 ./examples/raw 文件夹中，其中的文件如下所示：

.
|-- 2qbr.pdb
`-- 2qbr_ligand.sdf0 directories, 2 files

Interformer 是一个蛋白-配体复合物结构预测模型。第一个应用是为每个蛋白-配体原子对预测其相互作用感知的能量函数。这些能量函数可以用在传统的蒙特卡洛对接采样方法中，以生成高质量和合理的结合构象。

该模型的第二个应用涉及使用基于对比学习的姿势敏感亲和力和姿势评分模块。该模块为生成的对接姿势分配置信度分数，并预测相应的亲和力值。值得注意的是，对接姿势的质量越低，预测到的亲和力值（pIC50）越低的可能性也越大。

根据说明文档的步骤对配体和蛋白质数据进行准备。使用 OpenBabel 向配体添加氢原子，并确定其质子化状态。接下来，使用 RDKit 的 UFF 力场生成初始配体构象。下面是使用 Interformer 的案例测试。

激活运行环境，创建 ligand 文件夹以存放配体处理的文件。

conda activate interformer
mkdir -p examples/ligand

小分子准备

使用 openbabel 对配体质子化，添加氢原子。处理好的配体保存为 examples/ligand/2qbr_docked.sdf ，命令如下：

obabel examples/raw/2qbr_ligand.sdf \-p 7.4 \-O examples/ligand/2qbr_docked.sdf

创建 uff 文件夹，并使用 UFF 力场生成初始配体构象

mkdir -p examples/uffpython tools/rdkit_ETKDG_3d_gen.py \examples/ligand/ \examples/uff

UFF 力场生成初始配体构象保存为 examples/uff/2qbr_uff.sdf，运行输出如下，使用 UFF 力场生成初始配体构象和晶体结构（输入）的 RMSD 约为 1.87 Å 。

100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 9597.95it/s]
2qbr, RMSD:1.8653325928693114, E:76.88763102070008

蛋白&口袋准备

对于蛋白质，使用 Reduce 程序进行预处理，包括添加氢原子、质子化。然后，基于参考配体周围 10Å 半径的范围，得到蛋白口袋结构。

创建 pocket 文件夹，获取蛋白口袋输出的文件保存在该文件夹中。

mkdir -p examples/pocket

对整个蛋白质进行准备，包括加氢和质子化，输出结果保存为 examples/pocket/2qbr_reduce.pdb。运行命令如下：

reduce examples/raw/2qbr.pdb > examples/pocket/2qbr_reduce.pdb

运行输出如下：

Processing file: "examples/raw/2qbr.pdb"
Building or keeping OH & SH Hydrogens.
ERROR CTab(/usr/local/reduce_wwPDB_het_dict.txt): could not open
*WARNING*: Res "910" not in HETATM Connection Database. Hydrogens not added.
VDW dot density = 16/A^2
Probe radius = 0.25A
Orientation penalty scale = 1 (100%)
Eliminate contacts within 3 bonds.
Ignore atoms with |occupancy| <= 0.01 during adjustments.
Waters ignored if B-Factor >= 40 or |occupancy| < 0.66
Aromatic rings in amino acids accept hydrogen bonds.
Rotating NH3 Hydrogens.
Not processing Met methyls.Singles(size 62): A   5 LYS NZ  : A  12 LYS NZ  : A  20 TYR OH  : A  28 SER OG 
...
Adjusted 74 group(s)
If you publish work which uses reduce, please cite:
Word, et. al. (1999) J. Mol. Biol. 285, 1735-1747.
For more information see http://kinemage.biochem.duke.edu

使用项目提供提取蛋白口袋的脚本（tools/extract_pocket_by_ligand.py），需要传入三个参数，分别是 reduce 处理好的蛋白结构所在的文件夹路径（[protein_root]）、openbabel 处理好的配体结构所在的文件夹路径（[ligand_root]）以及是否从蛋白结构中删除参考配体结构，0 和 1 分别表示不删除和删除参考配体。

指定蛋白和配体结构文件夹分别为 examples/pocket/ 和 examples/ligand/，从蛋白结构中删除参考配体，具体命令如下：

python tools/extract_pocket_by_ligand.py \examples/pocket/ \examples/ligand/ \1

命令成功执行，输出如下，口袋结构保存到 examples/pocket/output/2qbr_pocket.pdb 。

[protein_root] [ligand_root] [0|1, rm_ccd or not]
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 92.25it/s]

为方便后续处理，把 examples/pocket/output/2qbr_pocket.pdb 移动到 examples/pocket 文件夹中。命令如下：

mv examples/pocket/output/2qbr_pocket.pdb examples/pocket

注：上述对原始配体和蛋白结构处理的步骤也可以使用其他 CADD 软件完成，模型的预测质量和处理步骤有直接关系。因此，任何在数据预处理过程中的提升都有可能增强模型的表现性能。

能量预测

准备好的文件保存在 ./example 文件夹中，模型输入需要特定的文件夹结构，目录结构如下：

.
|-- demo_dock.csv
|-- ligand
|   `-- 2qbr_docked.sdf
|-- pocket
|   |-- 2qbr_pocket.pdb
|   `-- 2qbr_reduce.pdb
|-- raw
|   |-- 2qbr.pdb
|   `-- 2qbr_ligand.sdf
`-- uff`-- 2qbr_uff.sdf4 directories, 7 files

其中，2qbr_docked.sdf 是晶体结构质子化并添加氢原子处理后的结构。2qbr_reduce.pdb 是通过 reduce 程序预处理后的整个蛋白结构，2qbr_pocket.pdb 则是根据配体提取出的小分子 7Å 距离范围的口袋结构。raw 文件夹中是原始的配体和蛋白结构。2qbr_uff.sdf 是 UFF 力场生成的初始配体构象。

项目提供 ./examples/demo_dock.csv 表格，程序将通过 Target 列中的 PDB ID 定位 ligand/、pocket/ 和 uff/ 文件夹中的文件，Molecule ID 列是配体的名字，该表格内容如下（注：如果没有这个 csv 文件则需要新建）：

数据预处理后，执行以下命令, 调用 Interfomer 模型生成预测能量函数文件：

DOCK_FOLDER=energy_outputPYTHONPATH=interformer/ python inference.py -test_csv examples/demo_dock.csv \
-work_path examples/ \
-ensemble checkpoints/v0.2_energy_model \
-batch_size 1 \
-posfix *val_loss* \
-energy_output_folder $DOCK_FOLDER \
-reload \
-debug

-work_path 指定输入数据的文件夹。-ensemble 指定训练好的模型路径。-batch_size 设置批次大小。-energy_output_folder 指定输出结果保存路径。

预测命令打印内容如下：

$$$$Interformer_Gnina2$$$$
Available GPUS:1
# Using GPUS:1
{'model': 'Interformer', 'debug': True, 'data_path': '', 'work_path': 'examples/', 'split_folder': 'diffdock_splits', 'ligand_folder': 'ligand', 'pocket_path': 'pocket', 'use_mid': False, 'affinity_pre': False, 'batch_size': 1, 'worker': 8, 'n_jobs': 70, 'method': 'Gnina2', 'Code': 'Flower', 'dataset': 'PLI', 'seed': 8888, 'clip': 5.0, 'checkpoint': './', 'gpus': 1, 'patience': 10, 'early_stop_metric': 'val_loss', 'early_stop_mode': 'min', 'num_epochs': 2000, 'main_loop': 1, 'precision': 16, 'per_target_sampler': 0, 'target_wise_sampler': 0, 'native_sampler': 0, 'filter_type': 'normal', 'use_ff_ligands': 'uff', 'num_ff_mols': 1, 'neg_weight': 1.0, 'num_nodes': 1, 'num_gpus': 8, 'per_target': False, 'metric_name': ['r'], 'test_csv': ['examples/demo_dock.csv'], 'reload': False, 'ensemble': ['checkpoints/v0.2_energy_model'], 'inference': False, 'posfix': '*val_loss*', 'energy_output_folder': 'energy_output', 'uff_as_ligand': False, 'n_layers': 6, 'num_heads': 8, 'hidden_dim': 128, 'ffn_scale': 4, 'intput_dropout_rate': 0.0, 'dropout_rate': 0.2, 'attention_dropout_rate': 0.2, 'attention_kd': 0, 'rbf_K': 128, 'rbf_cutoff': 10.0, 'intra_encoder': True, 'inter_encoder': True, 'attention_mask': True, 'edge_output_layer': True, 'weight_decay': 0.0001, 'warmup_updates': 1000, 'peak_lr': 0.0004, 'end_lr': 1e-09, 'pose_sel_mode': False, 'energy_mode': False, 'jizhi': False, 'exp': 'Interformer_Gnina2', 'node_featurizer': <feats.gnina_types.gnina_featurizer.PLIPAtomFeaturizer object at 0x76d69daf0e60>, 'edge_featurizer': <feats.gnina_types.gnina_featurizer.PLIPEdgeFeaturizer object at 0x76d69a8e7800>, 'complex_to_data': <function obabel_mol_parser at 0x76d69aaffb00>, 'angle_feat_size': 0, 'node_feat_size': 1, 'edge_feat_size': 1}
****************************************************************************************************
# Using Model <- checkpoints/v0.2_energy_model/checkpoints/general_PL_2020-epoch=085-val_loss=0.0191.ckpt
# [Eval] Using Affinity Reg Task.
# [Interformer] energy_output_fodler:
"args": {'model': 'Interformer', 'debug': False, 'data_path': '/apdcephfs_qy3/share_2932069/revoli/inter/csv/general_PL_2020.csv', 'work_path': '/apdcephfs_qy3/share_2932069/revoli/inter', 'split_folder': 'diffdock_splits', 'ligand_folder': 'ligand', 'pocket_path': 'pocket', 'use_mid': False, 'batch_size': 3, 'worker': 5, 'n_jobs': 70, 'method': 'Gnina2', 'Code': 'EnergyFix', 'dataset': 'sbdd', 'seed': 7461, 'clip': 5.0, 'checkpoint': '/apdcephfs_qy3/share_2932069/revoli/inter/beta', 'gpus': 8, 'patience': 20, 'early_stop_metric': 'val_loss', 'early_stop_mode': 'min', 'num_epochs': 2000, 'main_loop': 1, 'precision': 16, 'per_target_sampler': 0, 'target_wise_sampler': 0, 'native_sampler': 0, 'filter_type': 'normal', 'use_ff_ligands': 'uff', 'num_ff_mols': 1, 'mix_loss': False, 'neg_weight': 1.0, 'num_nodes': 1, 'num_gpus': 8, 'per_target': False, 'metric_name': ['r'], 'test_csv': ['./test.csv'], 'reload': True, 'ensemble': ['checkpoint'], 'inference': False, 'posfix': '*val_loss*', 'all_poses': False, 'energy_output_folder': '', 'uff_as_ligand': False, 'n_layers': 6, 'num_heads': 8, 'hidden_dim': 128, 'ffn_scale': 4, 'intput_dropout_rate': 0.0, 'dropout_rate': 0.1, 'attention_dropout_rate': 0.1, 'attention_kd': 0, 'rbf_K': 128, 'rbf_cutoff': 10.0, 'weight_decay': 1e-05, 'warmup_updates': 11000, 'peak_lr': 0.0012, 'end_lr': 1e-09, 'pose_sel_mode': False, 'energy_mode': True, 'jizhi': True, 'exp': 'Interformer_Gnina2', 'node_featurizer': <feats.gnina_types.gnina_featurizer.PLIPAtomFeaturizer object at 0x76d69a7128d0>, 'edge_featurizer': <feats.gnina_types.gnina_featurizer.PLIPEdgeFeaturizer object at 0x76d69a729700>, 'complex_to_data': <function obabel_mol_parser at 0x76d69aaffb00>, 'angle_feat_size': 0, 'node_feat_size': 1, 'edge_feat_size': 1}
$$$$Interformer_Gnina2$$$$
Number of parameters = 3,621,726
/workspace/anaconda3/envs/interformer/lib/python3.12/site-packages/lightning_fabric/connector.py:563: `precision=16` is supported for historical reasons but its usage is discouraged. Please set your precision to 16-mixed instead!
Using 16bit Automatic Mixed Precision (AMP)
GPU available: True (cuda), used: True
TPU available: False, using: 0 TPU cores
IPU available: False, using: 0 IPUs
HPU available: False, using: 0 HPUs
...
LOCAL_RANK: 0 - CUDA_VISIBLE_DEVICES: [0]
/workspace/anaconda3/envs/interformer/lib/python3.12/site-packages/pytorch_lightning/trainer/connectors/data_connector.py:441: The 'predict_dataloader' does not have many workers which may be a bottleneck. Consider increasing the value of the `num_workers` argument` to `num_workers=15` in the `DataLoader` to improve performance.
Predicting DataLoader 0:   0%|                                                                                                                                                                                                                                                                  | 0/1 [00:00<?, ?it/s][Interformer] Gussian Score ->2qbr-0
Predicting DataLoader 0: 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00,  3.65it/s]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
====================================================================================================
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
# Selected Model index..
/workspace/anaconda3/envs/interformer/lib/python3.12/site-packages/numpy/core/fromnumeric.py:3504: RuntimeWarning: Mean of empty slice.return _methods._mean(a, axis=axis, dtype=dtype,
/workspace/anaconda3/envs/interformer/lib/python3.12/site-packages/numpy/core/_methods.py:129: RuntimeWarning: invalid value encountered in scalar divideret = ret.dtype.type(ret / rcount)
checkpoints/v0.2_energy_model -> nan -> 
====================================================================================================
# Calculating Ensemble results....

（运行过程中似乎有报错，作者在项目的 issue 回答中提到，该报错是亲和力预测过程中计算 Pearson 相关分数时的报错，我们使用能量模型生成对接构象，可以忽略这个信息。）

模型预测速度很快，几秒钟就生成预测的能量函数。batch_size 设置为 1 ，GPU 最大占用约为 1 GB 。结果保存在 ./energy_output 文件夹中，生成文件如下：

.
|-- complex
|   `-- 2qbr_complex.pdb
|-- gaussian_predict
|   |-- 2qbr_G.db.bak
|   |-- 2qbr_G.db.dat
|   `-- 2qbr_G.db.dir
|-- ligand
|   `-- 2qbr_docked.sdf
`-- uff`-- 2qbr_uff.sdf4 directories, 6 files

其中，complex 中是通过参考配体截取的口袋。gaussian_predict 中是预测的能量函数文件。ligand 中的结构是从 examples/ 的 ligand 文件夹复制过来的，用于定位参考配体的对接盒子（20Å x 20Å x 20Å 采样空间）。uff 中的结构是从 examples/ 的 uff 文件夹复制过来的。

对接 pose 采样

接着，通过蒙特卡洛采样给定能量文件生成（采样）对接构象。设置 CPU 和 GPU 计算线程数目为 64，基于能量函数文件重构配体结构。尝试查找所有的对接构象，具体命令如下：

OMP_NUM_THREADS="64,64" python docking/reconstruct_ligands.py \-y \--cwd $DOCK_FOLDER \-y \--find_all find

运行输出如下：

 :: before initialization :: generators_for_threads.size() => [ 0 ] :: :: threadlimit = [ 2147483647 ] :: :: threadmax   = [ 64 ] :: :: i => [ 0 ] :: i_seed to i_generator_for_thread (pseudo-random generator) => [ 2098255030 ] :: :: i => [ 1 ] :: i_seed to i_generator_for_thread (pseudo-random generator) => [ 1634298443 ] :: :: i => [ 2 ] :: i_seed to i_generator_for_thread (pseudo-random generator) => [ 2449949960 ] :: 
...filename_input pdb_id  pose_rank  ...  intra_energy      rmsd sdf_rank
0   2qbr_docked.sdf   2qbr          0  ...     27.832400   1.12612        0
1   2qbr_docked.sdf   2qbr          1  ...     32.918147   5.40413        0
...
19  2qbr_docked.sdf   2qbr         19  ...     16.955457   9.25992        0
20  2qbr_docked.sdf   2qbr         20  ...     14.220284   0.00000        0[21 rows x 11 columns]

运行结束后，在 .//energy_output/ligand_reconstructing 路径生成了 20 个 pose （2qbr_docked.sdf 文件），文件夹内文件包括：

.
|-- 2qbr_docked.sdf
`-- 2qbr_docked.sdf_stat.csv0 directories, 2 files

2qbr_docked.sdf_stat.csv 文件则包含了，生成构象的能量、打分、排名信息，如下：

接着统计重构的对接构象的生成质量，默认生成构象和原始配体结构进行比较，如果使用 UFF 配体作为初始配体构象，可以添加此参数 --uff_folder uff

python docking/reconstruct_ligands.py \--cwd $DOCK_FOLDER \--find_all stat

统计数据保存为 ./energy_output/stat_ligand_reconstructing.csv，主要包括 rmsd、能量、torsions 数量和 pose_rank 等信息，表格部分内容如下（内容和 2qbr_docked.sdf_stat.csv 一致）：

把对接构象的统计信息表格和输入信息表格合并到一个表格中，以方便查看汇总数据。并把将对接后的 sdf 移动到工作路径（examples/infer，注：这里仅仅是结果合并，没有产生新的结果）。

python docking/merge_summary_input.py \$DOCK_FOLDER/ligand_reconstructing/stat_concated.csv \examples/demo_dock.csvcp -r $DOCK_FOLDER/ligand_reconstructing/ examples/infer

根据能量文件采样的分子对接构象以及排名、能量、RMSD 等信息的表格保存在 ./examples/infer 中，如下：

.
|-- 2qbr_docked.sdf
|-- 2qbr_docked.sdf_stat.csv
`-- stat_concated.csv0 directories, 3 files

预测 PoseScore & 亲和力

根据生成的分子构象，预测其 PoseScore 和亲和力值。

PYTHONPATH=interformer/ python inference.py \
-test_csv examples/demo_dock.round0.csv \
-work_path examples/ \
-ligand_folder infer/ \
-ensemble checkpoints/v0.2_affinity_model/model* \
-gpus 1 \
-batch_size 20 \
-posfix *val_loss* \
--pose_sel True

-work_path 指定预处理数据的路径。-ligand_folder 指定生成的对接构象的路径。-ensemble 指定使用的训练好的模型的路径。-batch_size 设置为 20 。

命令运行输出：

$$$$Interformer_Gnina2$$$$
Available GPUS:1
# Using GPUS:1
{'model': 'Interformer', 'debug': False, 'data_path': '', 'work_path': 'examples/', 'split_folder': 'diffdock_splits', 'ligand_folder': 'infer/', 'pocket_path': 'pocket', 'use_mid': False, 'affinity_pre': False, 'batch_size': 20, 'worker': 8, 'n_jobs': 70, 'method': 'Gnina2', 'Code': 'Flower', 'dataset': 'PLI', 'seed': 8888, 'clip': 5.0, 'checkpoint': './', 'gpus': 1, 'patience': 10, 'early_stop_metric': 'val_loss', 'early_stop_mode': 'min', 'num_epochs': 5, 'main_loop': 1, 'precision': 16, 'per_target_sampler': 0, 'target_wise_sampler': 0, 'native_sampler': 0, 'filter_type': 'normal', 'use_ff_ligands': 'uff', 'num_ff_mols': 1, 'neg_weight': 1.0, 'num_nodes': 1, 'num_gpus': 8, 'per_target': False, 'metric_name': ['r'], 'test_csv': ['examples/demo_dock.round0.csv'], 'reload': True, 'ensemble': ['checkpoints/v0.2_affinity_model/model0', 'checkpoints/v0.2_affinity_model/model1', 'checkpoints/v0.2_affinity_model/model2', 'checkpoints/v0.2_affinity_model/model3'], 'inference': False, 'posfix': '*val_loss*', 'energy_output_folder': '', 'uff_as_ligand': False, 'n_layers': 6, 'num_heads': 8, 'hidden_dim': 128, 'ffn_scale': 4, 'intput_dropout_rate': 0.0, 'dropout_rate': 0.2, 'attention_dropout_rate': 0.2, 'attention_kd': 0, 'rbf_K': 128, 'rbf_cutoff': 10.0, 'intra_encoder': True, 'inter_encoder': True, 'attention_mask': True, 'edge_output_layer': True, 'weight_decay': 0.0001, 'warmup_updates': 1000, 'peak_lr': 0.0004, 'end_lr': 1e-09, 'pose_sel_mode': True, 'energy_mode': False, 'jizhi': False, 'exp': 'Interformer_Gnina2', 'node_featurizer': <feats.gnina_types.gnina_featurizer.PLIPAtomFeaturizer object at 0x712870189a60>, 'edge_featurizer': <feats.gnina_types.gnina_featurizer.PLIPEdgeFeaturizer object at 0x712870189fd0>, 'complex_to_data': <function obabel_mol_parser at 0x712870337920>, 'angle_feat_size': 0, 'node_feat_size': 1, 'edge_feat_size': 1}
****************************************************************************************************
# Using Model <- checkpoints/v0.2_affinity_model/model0/checkpoints/general_PL_2020_round0_full-epoch=018-val_loss=0.4251.ckpt
# [Eval] Using Affinity Reg Task.
# [Eval] Using Pose Selection CLS task.
"args": {'model': 'Interformer', 'debug': False, 'data_path': '/apdcephfs_qy3/share_2932069/revoli/inter/csv/general_PL_2020_round0_full.csv', 'work_path': '/apdcephfs_qy3/share_2932069/revoli/inter', 'split_folder': 'diffdock_splits', 'ligand_folder': 'ligand', 'pocket_path': 'pocket', 'use_mid': False, 'affinity_pre': True, 'batch_size': 10, 'worker': 5, 'n_jobs': 70, 'method': 'Gnina2', 'Code': 'Affinity_230', 'dataset': 'sbdd', 'seed': 8741, 'clip': 5.0, 'checkpoint': '/apdcephfs_qy3/share_2932069/revoli/inter/beta', 'gpus': 8, 'patience': 20, 'early_stop_metric': 'val_pose_selection', 'early_stop_mode': 'max', 'num_epochs': 2000, 'main_loop': 1, 'precision': 16, 'per_target_sampler': 0, 'target_wise_sampler': 0, 'native_sampler': 0, 'filter_type': 'full', 'use_ff_ligands': 'uff', 'num_ff_mols': 1, 'neg_weight': 1.0, 'num_nodes': 1, 'num_gpus': 8, 'per_target': False, 'metric_name': ['r'], 'test_csv': ['./test.csv'], 'reload': True, 'ensemble': ['checkpoint'], 'inference': False, 'posfix': '*val_loss*', 'energy_output_folder': '', 'uff_as_ligand': False, 'n_layers': 6, 'num_heads': 8, 'hidden_dim': 128, 'ffn_scale': 4, 'intput_dropout_rate': 0.0, 'dropout_rate': 0.1, 'attention_dropout_rate': 0.1, 'attention_kd': 0, 'rbf_K': 128, 'rbf_cutoff': 10.0, 'intra_encoder': True, 'inter_encoder': True, 'attention_mask': True, 'edge_output_layer': True, 'weight_decay': 1e-05, 'warmup_updates': 10000, 'peak_lr': 0.0003, 'end_lr': 1e-09, 'pose_sel_mode': True, 'energy_mode': False, 'jizhi': True, 'exp': 'Interformer_Gnina2', 'node_featurizer': <feats.gnina_types.gnina_featurizer.PLIPAtomFeaturizer object at 0x71286ffd2b70>, 'edge_featurizer': <feats.gnina_types.gnina_featurizer.PLIPEdgeFeaturizer object at 0x71286c0dda90>, 'complex_to_data': <function obabel_mol_parser at 0x712870337920>, 'angle_feat_size': 0, 'node_feat_size': 1, 'edge_feat_size': 1}
$$$$Interformer_Gnina2$$$$
Number of parameters = 3,554,966
# csv <-- examples/demo_dock.round0.csv
[Bindingdata] Dataloader:lmdb
[Bindingdata] preprocess for affinity mode.
[Dataset] Loading previously saved cached file. <-- examples//tmp_beta/demo_dock.round0.csv-Gnina2-normal-uff--affinity-
@LMDB Loading <-examples//tmp_beta/demo_dock.round0.csv-Gnina2-normal-uff--affinity-
# [Lmdb] Geting Sampler Info <- examples//tmp_beta/demo_dock.round0.csv-Gnina2-normal-uff--affinity-.csv
# [Lmdb] Finished Loading Sampler Info.
[Bindingdata] Total Samples:21
...
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
# Selected Model index..
checkpoints/v0.2_affinity_model/model2 -> 1.000000 -> demo_dock:1.00
checkpoints/v0.2_affinity_model/model1 -> 1.000000 -> demo_dock:1.00
checkpoints/v0.2_affinity_model/model0 -> 1.000000 -> demo_dock:1.00
checkpoints/v0.2_affinity_model/model3 -> -1.000000 -> demo_dock:-1.00
====================================================================================================
# Calculating Ensemble results....
result/demo_dock.round0_ensemble.csv
Year             NaN
Resolution       NaN
pIC50            NaN
Molecule ID      NaN
pose_rank        NaN
num_torsions     NaN
energy           NaN
rmsd             NaN
pred_pIC50       NaN
pred_pose        NaN
pred_pIC50_var   NaN
Name: pIC50, dtype: float64----------------------------------------------------------------------------------------------------

命令运行，GPU 显存占用最大约为 10 GB，花费时间约 1 分钟。根据生成的分子构象，预测的 PoseScore 和亲和力值保存到 result/demo_dock.round0_ensemble.csv 表格中。

可以通过下面命令快速查看最终结果：

cat result/demo_dock.round0_ensemble.csv  # 可以查看最终结果

该表格部分内容如下，其中，energy 是每对蛋白质-配体原子的预测能量总和，pred_pIC50 是预测的亲和力值的负对数值，pred_pose 是预测的 PoseScore（对接构象的置信度评分）。该表格一共记录了这 21 个生成的分子构象情况。

Interformer 生成构象的排名前 3 的分子构象在蛋白口袋中的位置如下，蓝色的是晶体结构。

前三个对接构象和晶体结构的 RMSD、置信度打分、能量、预测的亲和力（pred_pIC50）和 PoseBusters 通过情况如下所示：

Rank	RMSD	置信度打分PoseScore	能量	pred_pIC50	PoseBusters check
0	0.85	0.97	-1348.26	5.91	True
1	5.18	0.81	-1230.33	5.18	True
2	3.66	0.55	-1204.99	3.66	True

可以看出，Interformer 模型判断好的构象并不只是参考 RMSD ，看重配体和蛋白之间的适应程度（能量和化学结构相容性）。

3.2.2 自定义案例（ABL1）

我们把 ABL1 体系作为自定义案例。6HD6.pdb 晶体结构，属于 Holo 结构，在此结构中，存在两个小分子，一个底物小分子抑制剂 Imatinib，另一个是变构小分子 Compound 6。除此以外，我们还找到 Compound 5 和 Compound 7 这两个变构分子，且活性与 Compound 6 有明显的区别。Compound N 则是来源于非激酶体系的变构小分子，是肯定不能结合在底物或者变构位点的。关于 ABL1 测试体系的详细介绍，可以参考前面介绍 FlexPose 的文章。分子的 2D 结构及活性情况如下所示：

Compound 6

Compound 6 排名前 3 （该排名和 PoseBuster 检查无关）的分子构象在蛋白口袋中的位置如下，蓝色的是晶体结构。可以看出在口袋内的片段和晶体结构相比变化不大，但接近溶剂区的片段摆动较大。

Compound 6 前三个对接构象和晶体结构的 RMSD、置信度打分、能量、预测的亲和力（pred_pIC50）和 PoseBusters 通过情况如下所示：

Rank	RMSD	置信度打分	能量	pred_pIC50	PoseBusters check
0	0.78	0.86	-958.29	5.91	True
1	0.74	0.58	-927.02	5.87	True
2	0.65	0.65	-901.66	3.64	True

Compound 5

Compound 5 排名前 3 （该排名和 PoseBuster 检查无关）的分子构象在蛋白口袋中的位置如下，蓝色的是晶体结构。可以看出前两个构象和晶体结构的位置基本一致，第三个构象的位置偏移较大。

Compound 5 前三个对接构象和晶体结构（cpd6）的最大公共子片段 RMSD、置信度打分、能量、预测的亲和力（pred_pIC50）和 PoseBusters 通过情况如下所示。可以看出这三个构象的置信度打分是比较低的。

Rank	RMSD	置信度打分	能量	pred_pIC50	PoseBusters check
0	0.65	0.33	-885.88	5.99	True
1	0.70	0.58	-851.64	6.41	True
2	1.00	0.06	-348.15	3.86	True

Compound 7

Compound 7 排名前 3 的分子构象在蛋白口袋中的位置如下，蓝色的是晶体结构。可以看出前三个构象和晶体结构的位置基本一致。

Compound 7 前三个对接构象和晶体结构（cpd6）的最大公共子片段 RMSD、置信度打分，能量、预测的亲和力（pred_pIC50）和 PoseBusters 通过情况如下所示。可以看出这三个构象的置信度打分是比较低的。

Rank	RMSD	置信度打分	能量	pred_pIC50	PoseBusters check
0	0.57	0.35	-945.65	6.82	True
1	0.57	0.02	-903.97	5.44	True
2	0.62	0.02	-849.41	5.49	True

Compound N

Compound N 排名前 3 的分子构象在蛋白口袋中的位置如下，蓝色的是晶体结构。可以看出前三个构象和晶体结构的位置和朝向基本一致。

Compound N 前三个对接构象和晶体结构（cpd6）的最大公共子片段 RMSD、置信度打分，能量、预测的亲和力（pred_pIC50）和 PoseBusters 通过情况如下所示。可以看出这三个构象的置信度打分是比较低的。

Rank	RMSD	置信度打分	能量	pred_pIC50	PoseBusters check
0	1.42	0.20	-709.95	5.30	True
1	1.45	0.13	-762.47	4.87	True
2	1.52	0.11	-602.44	4.77	True

下表汇总了自定义案例（ABL1）的四个分子 Top 1 的构象和晶体结构（cpd6）的最大公共子片段RMSD、置信度打分，能量、预测的亲和力（pred_pIC50）和 PoseBusters 通过情况如下所示。

化合物	RMSD	置信度打分	能量	IC50 (uM)	Kd (uM)	pred_pIC50	PoseBusters check
cpd5	0.65	0.33	-885.88	/	10	5.99(1.02)	True
cpd6	0.78	0.86	-958.29	550	2	5.91(1.23)	True
cpd7	0.57	0.35	-945.65	0.0023	<2	6.82(0.15)	True
cpdN	1.42	0.20	-709.95	/	/	5.30(5.01)	True

注：活性数值列中的 / 表示该化合物没有活性。pred_pIC50 列为模型预测的 IC50 的负对数值，括号中是换算成 IC50 的数值(单位是 uM ) ，以方便和实验活性比较。

3.3 训练 Interformer 测试

项目提供处理好的训练数据，保存在 zenodo，网址为：Interformer: An Interaction-Aware Model for Protein-Ligand Docking and Affinity Prediction 。找到对应的下载位置（如下所示），下载 interformer_train.tar.gz 压缩文件。

下载后上传到 ./ ，通过以下命令解压文件：

tar -xzvf interformer_train.tar.gz

解压得到 ./poses 和 ./train 两个文件夹，./poses 文件夹中保存的是结构文件，./train 文件夹中保存的是数据信息表格。

项目提供训练的脚本 ./scripts/train/run_single.sh ，我们的机器只有一个 GPU，指定 CUDA_VISIBLE_DEVICES=0，并设置 -gpus 为 1，batch_size 设置为 4，具体训练命令如下：

PYTHONPATH=interformer/ CUDA_VISIBLE_DEVICES=0 python train.py -data_path /workspace/GuanXL/projects/interformer/train/general_PL_2020.csv \
-work_path /workspace/GuanXL/projects/interformer/poses \
-dataset sbdd \
-seed 1111 \
-filter_type normal \
-per_target_sampler 0 \
-native_sampler 1 \
-Code Normal \
--warmup_updates \
10000 \
--peak_lr \
0.0006 \
-model \
Interformer \
-batch_size \
4 \
-gpus \
1 \
-method \
Gnina2 \
--n_layers \
6 \
--hidden_dim \
128 \
--num_heads \
8 \
--dropout_rate \
0.1 \
--attention_dropout_rate \
0.1 \
--weight_decay \
1e-5 \
-patience 30 \
-early_stop_metric val_loss \
-early_stop_mode min \
-affinity_pre

我们在后台运行该训练脚本，命令如下，输出记录在 train.log 中，以便查看。

nohup bash -c "time bash scripts/train/run_single.sh" > train.log 2>&1 &

我们测试训练了 5 个 epoch 之后，batch_size 设置为 4 ，GPU 显存最大占用约 20 GB 。输出默认保存在 ./lightning_logs/general_PL_2020_Interformer_Normal 文件夹中。训练好的模型文件保存在 ./lightning_logs/general_PL_2020_Interformer_Normal/version_3/checkpoints 文件夹中，如下所示：

.
|-- general_PL_2020-epoch=003-val_loss=1.1954.ckpt
|-- general_PL_2020-epoch=004-val_affinity=0.5325.ckpt
|-- last-v1.ckpt
`-- last.ckpt0 directories, 4 files

其中，general_PL_2020-epoch=003-val_loss=1.1954.ckpt 表示训练的这 5 个 epoch 中，按照验证集损失这个指标表现最好的模型是第 3 个 epoch，该文件是保存的模型文件。general_PL_2020-epoch=004-val_affinity=0.5325.ckpt 表示训练的这 5 个 epoch 中，按照验证集亲和力这个指标表现最好的模型是第 4 个 epoch，该文件是保存的模型文件。last.ckpt 是保存的最后一个 epoch 的模型文件。

四、总结

确定蛋白质-配体复合物的结构在药物开发领域中一直是一个重要的挑战。为了应对这一挑战，作者提出了 Interformer ，这是一种专为蛋白质-配体对接和亲和力预测设计的深度学习生成模型。该模型整合了一个强大的、具有相互作用感知能力的 MDF ，能够高效地恢复具体的相互作用。此外，Interformer 的机制易于解释，能够有效解决当前最先进的深度学习模型常常忽视重要非共价相互作用的常见不足。

除了阐明药物设计中结合模式的重要性外，根据配体对特定靶标的亲和力对其进行排序或筛选同样至关重要。考虑到许多最先进的深度学习模型容易在晶体结构上过拟合，Interformer 采用了一种训练策略，利用对比学习和负采样来增强对接构象的敏感性。这种方法使得 Interformer 能够通过关注蛋白质和配体原子对之间的具体相互作用，区分较不精确和更有利的对接构象。这一强大的能力使得该模型能够在实际场景中提高预测的泛化能力。

Interformer 在两个广泛认可的基准测试中展示了蛋白质-配体对接的持续改进，生成了物理上合理的对接构象，增强了下游应用的潜力。在亲和力预测领域，Interformer 在四个内部真实世界亲和力基准上也表现出持续的进步。进一步应用 Interformer 于两个内部药物开发管道，成功识别了两个具有纳摩尔级活性的小分子。