SQUID - 形状条件下的基于分子片段的3D分子生成等变模型评测

SQUID 是一个形状条件下基于片段的3D分子生成模型，给一个3D参考分子，SQUID 可以根据参考分子的形状，基于片段库，生成与参考分子形状非常相似的分子。

SQUID 模型来自于 ICLR 2023 文章（2022年10月6日提交），对应的arix文献为《Equivariant Shape-Conditioned Generation of 3D Molecules for Ligand-Based Drug Design》, 文章链接为：https://arxiv.org/abs/2210.04893

SQUID是: Shape-Conditioned Equivariant Generator for Drug-Like Molecules的缩写（神奇的脑回路），是麻省理工 Keir Adams 和 Connor W. Coley的工作，似乎两人研究领域都是AI4science。

一、模型介绍

1.1 基本介绍

基于形状的虚拟筛选在基于配体的药物设计中至关重要，主要目的是识别与已知配体具有相似 3D 形状的分子。传统方法（例如：openeye 的 ROCS）依赖于枚举的化学库，这限制了新化学空间的探索。

基于形状的 3D 分子生成模型，需要完成任意构象的成对比较、形状相似性度量、2D 图拓扑对 3D 形状的依赖以及可旋转键的灵活性处理等任务。为此，作者开发了SQUID模型。SQUID 是在形状条件下 3D 分子生成模型，可以用于在形状条件下的化学空间探索任务。

SQUID 模型将分子形状使用等变点云网络（ equivariant point cloud networks）编码成旋转不变形的的等变特征；使用图神经网路将2D的化学结构变分编码为化学的不变特征；然后，将化学的不变特征和形状的等变特征合并在一起，生成多样性的 3D 分子。

在 3D 分子的生成过程中，SQUID 使用基于片段的策略，固定局部键长和角度（人工分子），优先考虑可旋转键。大大简化了 3D 坐标生成，可以生成类药物分子，同时保持化学和几何有效性。SQUID 设计一个可旋转的键评分网络，学习局部键旋转如何影响整体形状，使我们的解码器能够生成最适合目标形状的 3D 构象。

SQUID 3D分子生成过程示意图：

SQUID 模型可以解开控制 2D 化学特性和 3D 形状的潜在变量，从而能够零样本生成与任何编码目标具有相似形状的拓扑不同分子。SQUID 模型既预测分子图，也进行分子的坐标预测，但是使用基于片段的策略，而且是在形状条件下，分子从无到有的生成。简单的来说，在片段库内，根据参考分子的3D形状，生成形状极为相似的分子。

在 github 中，作者给了一段视频，看起来模型非常有意思，可以在生成形状非常像的小分子：

SQUID_demo

1.2 模型框架介绍

SQUID 使用编码器-解码器架构，将 3D 分子形状和化学身份进行编码，并生成与目标形状相匹配的新分子。在分子图的表示方面，节点（原子）特征有：原子质量、原子编号的独热编码、原子电荷和芳香性的独热编码、单键、双键、芳香键和三键的数量独热编码；边的特征为键序的独热编码。在形状表示方面，通过 3D 分子中从每个原子中心的高斯分布中采样点云来表示分子的 3D 形状。

SQUID的编码器和解码器结构如下图：

特征的编码分为三块：（1）片段编码器：使用图神经网络（GNN）对片段进行编码，得到全局的片段嵌入。（2）形状编码器：使用 VN-DGCNN （Vector Neurons (VN)-based equivariant point cloud encoder VN-DGCNN）将分子的点云编码为等变的每点向量特征，并进行局部均值池化，得到每个原子的等变表示。（3）使用GNN将分子图编码为学习到的原子嵌入，同时条件化在不变的形状特征上。

在编码完成后，混合形状和变分特征，将变分原子特征通过线性变换与等变形状特征混合，并通过VN-MLP进行处理，生成全局的等变和不变表示。

在分子生成的过程（解码过程）中，从目标分子中提取一个小的子结构作为生成的起始结构。在生成过程中，每次新添加原子或片段之前，对部分分子进行编码，确保与目标形状对齐。在分子图上使用图生成器，预测是否停止生成、选择焦点原子、选择新片段、预测附着位点和键序。生成过程中，动作选择会遵循已知的化学价规则，以确保化学有效性。新生成的键使用键打分器，模型通过预测旋转角度ψfoc来最大化生成分子与目标形状的相似性。

因此，SQUID是一个基于片段的、形状条件下的、自回归的 3D 分子生成的 VAE 模型，目的是生成有效的、符合特定 3D 形状的分子。

关于 SQUID 模型的损失，训练模型损失的表达式如下：

$L_{total} = L_{graph-gen} + β_{KL} L_{KL} + β_{next-shape} L_{next-shape} + β_{∅-shape} L_{∅-shape} + L_{bond-scoring}$

其中， $L_{graph-gen}$ 代表分子图生成的损失，具体包括：分子可继续生长与否损失、生长点选择损失、下一个片段选择损失、下一个片段的附着位点损失、以及链接分子片段与下一个片段的键类型的损失，这些损失都是交叉熵损失。 $L_{KL}$ 代表散度损失，约束变分编码。 $L_{next-shape}$ 为辅助损失，用于引导生成器生成形状相似的多样化分子，表示生成分子与目标形状之间的匹配程度。 $L_{bond-scoring}$ 是旋转键的回归评分损失。β为各种损失的权重，是模型的超参数。

关于模型的数据集，模型使用来自MOSES数据集的药物样分子进行训练，最大原子数量为27。从数据集中提取了100个片段（涵盖24种原子类型），确保每个片段包含两个或两个以上的原子，主要为环状片段，如苯环、萘环等。然后以此为基础，生成包含 1.3M 个 3D 分子的最终数据集。

每个训练分子从其3-6个原子的 3D 子结构中提取一个子结构，作为生成的起始结构。对每个分子生成一个 3D 构象，设定非环键的键长为其经验均值，并使用启发式规则固定非环键角。虽然这种处理忽略了真实分子的键变形，但对全局形状影响较小，并且可以在不严重改变形状的情况下恢复精细几何结构。

1.3 模型性能介绍

论文对SQUID模型的性能进行了两个方面的评估，分别是：形状条件下的多样性分子生成和形状限制下的分子优化。

1.3.1 形状条件下的多样性分子生成

形状条件下的多样性分子生成的目的是：评估模型生成与目标形状相似但化学结构多样的分子的能力。采用的方法是，对于测试集中的1000个分子，使用不同的\lambda（分别为0.3和1.0，采样分布，代表后验知识的比例），每个分子生成50个候选分子。使用高斯描述计算的体积重叠，量化生成分子与目标分子的3D形状相似性（simS）；使用RDKit计算的Tanimoto相似性（simG），量化生成分子与目标分子的化学结构相似性。

生成的分子中。simG< 0.7 或者 0.3的分子的 simS 分布，结果如下图：

图中表明，使用λ=0.3时，生成的分子化学多样性较低，但形状相似性较高；使用λ=1.0时，生成的分子化学多样性较高，但形状相似性稍有下降。另外，其他统计表明，99%的生成分子是新颖的，95%是独特的，且100%化学有效，87%的生成分子没有任何空间冲突，表明SQUID生成的3D分子具有较高的真实性和化学有效性。从图C来看，生成分子的 3D 结构与目标3D结构之间非常相似。

1.3.2 形状限制下的分子优化

形状限制下的分子优化的目的是：SQUID 在保持高形状相似度的同时，优化分子的某些目标属性（如生物活性、毒性等）的能力。使用的方法是：使用GaucaMol基准中的目标属性，通过遗传算法对种子分子的变分编码进行迭代变异，生成目标分数更高且形状相似的分子。

SQUID 项目的的 github 连接： https://github.com/keiradams/SQUID

二、环境安装

复制项目代码：

git clone https://github.com/keiradams/SQUID.git

项目目录如下：

.
├── data # 下载的数据
├── dataset_generation
├── data.zip # 下载的数据
├── environment.yml
├── LICENSE
├── models
├── MO_virtual_screening
├── our_test
├── paper_results.zip # 下载的数据
├── README.md
├── RUN_ME.ipynb
├── shape_conditioned_generation_dataset_baseline.py
├── shape_conditioned_generation_evaluations.py
├── shape_constrained_optimization_evaluations.py
├── similarity-all_hits.txt
├── similarity-all.txt
├── test
├── trained_models
├── train_graph_generator.py
├── train_scorer.py
└── utils

创建conda环境：

conda env create -f environment.yml

在创建环境的过程中，会出现以下报错：

Pip subprocess error:
ERROR: Could not find a version that satisfies the requirement openeye-toolkits==2022.1.1 (from versions: none)
ERROR: No matching distribution found for openeye-toolkits==2022.1.1failedCondaEnvException: Pip failed

这是在安装 openeye-toolkits 找不到包的问题。openeye-toolkits 不是一个开源的免费工具，因此无法直接通过 pip 安装。

如果有 openeye-toolkits 安装包，可以通过 python setup.py install 安装（注：linse需要放置在安装包的 Other/Proprietary License 文件夹内）。执行完 python setup.py install ，pip 还会从 pipy 中安装其中的依赖，但是会找不到，因为openeye-toolkits 不是开源的，因此也找不到依赖，直接kill就好。

openeye 工具包使用 ROCS(Rapid Overlay of Chemi cal Structures) 将 3D 分子刚性对齐，并计算相似性。如果没有 openeye 的 linsece 可以使用 Shaep 代替，但是计算效率会下降。

Shaep 实现基于形状和静电势的分子结构刚性叠加功能。Shaep 的下载地址是：ShaEP。下载完成后，直接解压即可以用，无需安装。解压后放置在 ../目录中。

在使用 conda env create -f environment.yml 安装环境的时候，由于 openeye 包的安装失败，打断了其他模块的安装，会导致其他模块没有安装。因此需要继续安装：

conda install pyg pytorch-cluster pytorch-scatter -c pyg

三、SQUID 分子生成测试

3.1 Notebook中的示例

作者提供了一个 notebook （./RUN_ME.ipynb），在其中以示例，给出了使用方法介绍。

在运行之前，需要下载数据集。下载链接为：https://figshare.com/s/3d2f8fd57d9a65fe237e

下载后，解压到 ./ 目录中，会有一个data文件夹，比较大，共 47G。其中，包含模型训练和测试需要的 MOSES 数据集、100个分子片段组成的片段库、过滤后的分子库、RDKIT生成的小分子的3D构象以及固定烷基化学键距离和角度前后的构象、预训练图生成器和键打分器的数据。

解压完成后的 ./data文件目录为：

.
└── MOSES2├── max_future_rocs_data_artificial_alpha_2_0├── MOSES2_training_val_atom_fragment_associations_array.npy├── MOSES2_training_val_AtomFragment_database.pkl├── MOSES2_training_val_atoms_pointer.npy├── MOSES2_training_val_bond_lookup.pkl├── MOSES2_training_val_bonds_pointer.npy├── MOSES2_training_val_canonical_terminalSeeds_unmerged_future_rocs_database_all_reduced.pkl├── MOSES2_training_val_edge_features_array.npy├── MOSES2_training_val_edge_index_array.npy├── MOSES2_training_val_filtered_database_artificial_mols.pkl├── MOSES2_training_val_filtered_database.pkl├── MOSES2_training_val_node_features_array.npy├── MOSES2_training_val_unique_atoms.npy├── MOSES2_training_val_xyz_array.npy├── MOSES2_training_val_xyz_artificial_array.npy├── MOSES2_train_smiles_split.csv├── MOSES2_val_smiles_split.csv├── test_MOSES.csv├── test_MOSES_filtered_artificial_mols.pkl├── training_split_arrays├── train_MOSES.csv└── validation_split_arrays4 directories, 19 files

以下是 Notebook 中的内容。

导入模块：

import torch_geometric
import torch
import torch_scatter
import math
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from copy import deepcopy
import rdkit
import rdkit.Chem
import rdkit.Chem.AllChem
import rdkit.Chem.rdMolTransforms
import networkx as nx
import random
from tqdm import tqdm
from rdkit.Chem import rdMolTransforms
import itertools
import os
import pickleimport torch.nn as nn
import torch.nn.functional as F
from models.vnn.models.vn_layers import *
from models.vnn.models.utils.vn_dgcnn_util import get_graph_featurefrom utils.general_utils import *from models.EGNN import *
from models.models import *

模型的超参数，包括，先验插值比例 \lambda，局部生长停止的阈值、标准偏差程度，是否使用等变网络。

interpolate_to_GNN_prior = 1.0 # 'prior'
stop_threshold = 0.01
variational_GNN_factor = 1.0
ablateEqui = False

准备对10个分子进行测试，每个分子生成20个分子。

repetitions = 20
total_evaluations = 10

加载必要的数据，包括：片段库、键长表、原子种类字典、测试集的3D分子

# 使用人工的分子，即分子中烷基键的长固定
use_artificial_mols = TrueAtomFragment_database = pd.read_pickle('data/MOSES2/MOSES2_training_val_AtomFragment_database.pkl')
AtomFragment_database = AtomFragment_database.iloc[1:].reset_index(drop = True) # removing stop token from AtomFragment_database
fragment_library_atom_features = np.concatenate(AtomFragment_database['atom_features'], axis = 0).reshape((len(AtomFragment_database), -1))bond_lookup = pd.read_pickle('data/MOSES2/MOSES2_training_val_bond_lookup.pkl')
unique_atoms = np.load('data/MOSES2/MOSES2_training_val_unique_atoms.npy') test_mol_df = pd.read_pickle('data/MOSES2/test_MOSES_filtered_artificial_mols.pkl')

test_mol_df 是一个 dataframe 对象，如下图：

使用人工分子，然后打乱顺序

test_mols = list(test_mol_df.artificial_mol)random.seed(0)
indices_to_evaluate = list(range(0, len(test_mols)))
random.shuffle(indices_to_evaluate)val_mols = [test_mols[i] for i in indices_to_evaluate]
val_mols_index = indices_to_evaluate

图生成器的 checkpoint、键打分器的 checkpoint：

if not ablateEqui:model_3D_PATH = 'trained_models/graph_generator.pt'rocs_model_3D_PATH = 'trained_models/scorer.pt'
else:model_3D_PATH = 'trained_models/graph_generator_ablateEqui.pt'rocs_model_3D_PATH = 'trained_models/scorer_ablateEqui.pt'

初始化图生成器和键打分器，并加载权重：

# HYPERPARAMETERS for 3D graph generator
pointCloudVar = 1. / (12. * 1.7) model_3D = Model_Point_Cloud_Switched(input_nf = 45, edges_in_d = 5, n_knn = 5, conv_dims = [32, 32, 64, 128], num_components = 64, fragment_library_dim = 64, N_fragment_layers = 3, append_noise = False, N_members = 125 - 1, EGNN_layer_dim = 64, N_EGNN_layers = 3, output_MLP_hidden_dim = 64, pooling_MLP = False, shared_encoders = False, subtract_latent_space = True,variational = False,variational_mode = 'inv', # not usedvariational_GNN = True,mix_node_inv_to_equi = True,mix_shape_to_nodes = True,ablate_HvarCat = False,predict_pairwise_properties = False,predict_mol_property = False,ablateEqui = ablateEqui,old_EGNN = False,).float()model_3D.load_state_dict(torch.load(model_3D_PATH, map_location=next(model_3D.parameters()).device), strict = True)
model_3D.eval()# HYPERPARAMETERS for ROCS scorer
rocs_pointCloudVar = 1. / (12. * 1.7) rocs_model_3D = ROCS_Model_Point_Cloud(input_nf = 45, edges_in_d = 5, n_knn = 10, conv_dims = [32, 32, 64, 128], num_components = 64, fragment_library_dim = 64,N_fragment_layers = 3, append_noise = False, N_members = 125 - 1, EGNN_layer_dim = 64, N_EGNN_layers = 3, output_MLP_hidden_dim = 64, pooling_MLP = False, shared_encoders = False, subtract_latent_space = True,variational = False,variational_mode = 'inv', # not usedvariational_GNN = False,mix_node_inv_to_equi = True,mix_shape_to_nodes = True,ablate_HvarCat = False,ablateEqui = ablateEqui,old_EGNN = False,).float()rocs_model_3D.load_state_dict(torch.load(rocs_model_3D_PATH, map_location=next(rocs_model_3D.parameters()).device), strict = True)
rocs_model_3D.eval()

执行SQUID，生成分子

# muting noisy warnings
from rdkit import RDLogger
import warnings
lg = RDLogger.logger()
lg.setLevel(RDLogger.CRITICAL)
warnings.filterwarnings("ignore", category=DeprecationWarning) 
warnings.filterwarnings("ignore", category=UserWarning) 
warnings.filterwarnings("ignore", category=Warning) 
np.warnings.filterwarnings('ignore', category=np.VisibleDeprecationWarning)file_idx = 0reference_mols_list = []
unaligned_mols_list = []for m_idx, m_ in enumerate(val_mols):seed = 0random.seed(seed)np.random.seed(seed = seed)torch.manual_seed(seed)m = deepcopy(m_)mol = deepcopy(m)# 生成 3D构象 即坐标xyz = np.array(mol.GetConformer().GetPositions())# 中心center_of_mass = np.sum(xyz, axis = 0) / xyz.shape[0]# 去中心xyz_centered = xyz - center_of_mass# 计算每个原子的隐向量for i in range(0, mol.GetNumAtoms()):x,y,z = xyz_centered[i]mol.GetConformer().SetAtomPosition(i, Point3D(x,y,z))m = deepcopy(mol)mol_target = deepcopy(m)ring_fragments = get_ring_fragments(mol_target)all_possible_seeds = get_all_possible_seeds(mol_target, ring_fragments)terminal_seeds = filter_terminal_seeds(all_possible_seeds, mol_target)select_seeds = get_starting_seeds(mol_target, AtomFragment_database, fragment_library_atom_features, unique_atoms, bond_lookup)if len(select_seeds) == 0:continuerandom_seed_selection = random.randint(0, len(select_seeds) - 1)select_seeds = [select_seeds[random_seed_selection]] * repetitionsrepeated_rocs = []repeated_tanimoto = []reference_mols = []unaligned_mols = []for seed in select_seeds:mol = deepcopy(m)# extracting starting seed and preparing to generateframe_generation, frame_rocs = get_frame_terminalSeeds(mol, seed, AtomFragment_database, include_rocs = True)positions = list(frame_rocs.iloc[0].positions_before)start = 0for i in range(len(frame_generation)):if (set(frame_generation.iloc[i].partial_graph_indices) == set(positions)) & (frame_generation.iloc[i].next_atom_index == -1):start = i + 1breakif len(frame_generation.iloc[0].partial_graph_indices) == 1: # seed is a terminal ATOMterminalSeed_frame = frame_generation.iloc[0:start].reset_index(drop = True)sequence = get_ground_truth_generation_sequence(terminalSeed_frame, AtomFragment_database, fragment_library_atom_features)mol = deepcopy(terminalSeed_frame.iloc[0].rdkit_mol_cistrans_stereo)partial_indices = deepcopy(terminalSeed_frame.iloc[0].partial_graph_indices_sorted)final_partial_indices = deepcopy(terminalSeed_frame.iloc[-1].partial_graph_indices_sorted)ring_fragments = get_ring_fragments(mol)add_to_partial = [list(f) for p in final_partial_indices for f in ring_fragments if p in f]add_to_partial = [item for sublist in add_to_partial for item in sublist]final_partial_indices = list(set(final_partial_indices).union(add_to_partial))queue_indices = deepcopy(terminalSeed_frame.iloc[0].focal_indices_sorted)_, seed_mol, queue, positioned_atoms_indices, atom_to_library_ID_map, _, _, _ = generate_seed_from_sequence(sequence, mol, partial_indices, queue_indices, AtomFragment_database, unique_atoms, bond_lookup, stop_after_sequence = True)seed_node_features = getNodeFeatures(seed_mol.GetAtoms())for k in atom_to_library_ID_map:seed_node_features[k] = AtomFragment_database.iloc[atom_to_library_ID_map[k]].atom_featuresG = get_substructure_graph(mol, final_partial_indices)G_seed = get_substructure_graph(seed_mol, list(range(0, seed_mol.GetNumAtoms())), node_features = seed_node_features)nm = nx.algorithms.isomorphism.generic_node_match(['atom_features'], [None], [np.allclose])em = nx.algorithms.isomorphism.numerical_edge_match("bond_type", 1.0)GM = nx.algorithms.isomorphism.GraphMatcher(G, G_seed, node_match = nm, edge_match = em)assert GM.is_isomorphic()idx_map = GM.mappingelse: # seed is a terminal FRAGMENTpartial_indices = deepcopy(frame_generation.iloc[0].partial_graph_indices_sorted)final_partial_indices = partial_indicesseed_mol = generate_conformer(get_fragment_smiles(mol, partial_indices))idx_map = get_reindexing_map(mol, partial_indices, seed_mol)positioned_atoms_indices = sorted([idx_map[f] for f in final_partial_indices])atom_to_library_ID_map = {} # no individual atoms yet generatedqueue = [0] # 0 can be considered the focal root nodefor i in final_partial_indices:x,y,z = mol.GetConformer().GetPositions()[i]seed_mol.GetConformer().SetAtomPosition(idx_map[i], Point3D(x,y,z)) starting_queue = deepcopy(queue)try:_, updated_mol, _, _, _, N_rocs_decisions, _, _, _, chirality_scored = generate_3D_mol_from_sequence(sequence = [], partial_mol = deepcopy(seed_mol), mol = deepcopy(mol_target), positioned_atoms_indices = deepcopy(positioned_atoms_indices), queue = starting_queue, atom_to_library_ID_map = deepcopy(atom_to_library_ID_map), model = model_3D, rocs_model = rocs_model_3D,AtomFragment_database = AtomFragment_database,unique_atoms = unique_atoms, bond_lookup = bond_lookup,N_points = 5, N_points_rocs = 5,stop_after_sequence = False,mask_first_stop = False,stochastic = False, chirality_scoring = True,stop_threshold = stop_threshold,steric_mask = True,variational_factor_equi = 0.0,variational_factor_inv = 0.0, interpolate_to_prior_equi = 0.0,interpolate_to_prior_inv = 0.0, use_variational_GNN = True, variational_GNN_factor = variational_GNN_factor, interpolate_to_GNN_prior = interpolate_to_GNN_prior, rocs_use_variational_GNN = False, rocs_variational_GNN_factor = 0.0, rocs_interpolate_to_GNN_prior = 0.0,pointCloudVar = pointCloudVar, rocs_pointCloudVar = rocs_pointCloudVar,)pred_rocs = get_ROCS(torch.tensor(np.array(updated_mol.GetConformer().GetPositions())), torch.tensor(np.array(mol.GetConformer().GetPositions())))tanimoto = rdkit.DataStructs.FingerprintSimilarity(*[rdkit.Chem.RDKFingerprint(x) for x in [mol, updated_mol]])reference_mols.append(mol)unaligned_mols.append(updated_mol)repeated_rocs.append(pred_rocs.item())repeated_tanimoto.append(tanimoto)except Exception as e:print(f'error during 3D generation -- {e}')continueprint(f'Encoded Molecule {file_idx + 1}:')print('(non-aligned) shape scores:', [np.round(r, 3) for r in repeated_rocs])print('tanimoto chemical similarity:', [np.round(r, 3) for r in repeated_tanimoto])print()file_idx += 1reference_mols_list.append(reference_mols[0])unaligned_mols_list.append(unaligned_mols)if file_idx == total_evaluations:break

运行输出。生成的分子保存在 unaligned_mols_list，参考分子保存在 reference_mols_list。每一个分子没有经过 align 之前的形状相似性、化学结构相似性：

Encoded Molecule 1:
(non-aligned) shape scores: [0.45, 0.645, 0.55, 0.627, 0.617, 0.512, 0.435, 0.802, 0.37, 0.532, 0.438, 0.527, 0.778, 0.61, 0.553, 0.579, 0.822, 0.517, 0.807, 0.476]
tanimoto chemical similarity: [0.266, 0.253, 0.271, 0.352, 0.208, 0.232, 0.332, 0.465, 0.304, 0.257, 0.27, 0.372, 0.375, 0.244, 0.252, 0.289, 0.327, 0.348, 0.276, 0.255]
warning: bond distance between atoms 8 and 8 unknown
Encoded Molecule 2:
(non-aligned) shape scores: [0.421, 0.529, 0.492, 0.708, 0.486, 0.663, 0.68, 0.627, 0.693, 0.707, 0.582, 0.721, 0.784, 0.655, 0.633, 0.674, 0.513, 0.516, 0.704, 0.494]
tanimoto chemical similarity: [0.338, 0.416, 0.188, 0.316, 0.315, 0.215, 0.342, 0.27, 0.339, 0.291, 0.277, 0.302, 0.296, 0.285, 0.46, 0.278, 0.318, 0.266, 0.248, 0.195]
Encoded Molecule 3:
(non-aligned) shape scores: [0.431, 0.399, 0.382, 0.625, 0.512, 0.663, 0.414, 0.627, 0.748, 0.467, 0.568, 0.434, 0.686, 0.815, 0.573, 0.448, 0.571, 0.472, 0.576, 0.389]
tanimoto chemical similarity: [0.287, 0.172, 0.135, 0.254, 0.29, 0.207, 0.179, 0.252, 0.297, 0.142, 0.271, 0.16, 0.187, 0.265, 0.298, 0.274, 0.16, 0.306, 0.279, 0.265]
Encoded Molecule 4:
(non-aligned) shape scores: [0.251, 0.374, 0.785, 0.333, 0.371, 0.416, 0.533, 0.384, 0.365, 0.379, 0.361, 0.336, 0.499, 0.352, 0.356, 0.349, 0.472, 0.329, 0.602, 0.411]
tanimoto chemical similarity: [0.243, 0.24, 0.246, 0.2, 0.232, 0.223, 0.181, 0.206, 0.237, 0.223, 0.195, 0.161, 0.207, 0.227, 0.2, 0.194, 0.171, 0.162, 0.258, 0.179]
Encoded Molecule 5:
(non-aligned) shape scores: [0.439, 0.475, 0.649, 0.666, 0.328, 0.605, 0.574, 0.419, 0.345, 0.516, 0.701, 0.209, 0.233, 0.57, 0.666, 0.515, 0.528, 0.54, 0.631, 0.59]
tanimoto chemical similarity: [0.164, 0.147, 0.165, 0.172, 0.136, 0.218, 0.159, 0.167, 0.119, 0.157, 0.173, 0.18, 0.137, 0.169, 0.15, 0.148, 0.154, 0.159, 0.171, 0.152]
Encoded Molecule 6:
(non-aligned) shape scores: [0.72, 0.534, 0.828, 0.375, 0.767, 0.729, 0.818, 0.675, 0.65, 0.553, 0.607, 0.587, 0.767, 0.678, 0.444, 0.764, 0.753, 0.663, 0.625, 0.598]
tanimoto chemical similarity: [0.215, 0.236, 0.2, 0.169, 0.198, 0.181, 0.15, 0.21, 0.169, 0.232, 0.19, 0.185, 0.245, 0.164, 0.209, 0.215, 0.181, 0.166, 0.181, 0.208]
Encoded Molecule 7:
(non-aligned) shape scores: [0.367, 0.416, 0.445, 0.328, 0.443, 0.482, 0.397, 0.27, 0.198, 0.583, 0.464, 0.513, 0.52, 0.579, 0.53, 0.408, 0.418, 0.472, 0.239, 0.152]
tanimoto chemical similarity: [0.199, 0.238, 0.238, 0.197, 0.245, 0.248, 0.177, 0.246, 0.119, 0.2, 0.213, 0.227, 0.187, 0.283, 0.216, 0.222, 0.189, 0.322, 0.191, 0.227]
Encoded Molecule 8:
(non-aligned) shape scores: [0.525, 0.801, 0.493, 0.626, 0.485, 0.72, 0.585, 0.618, 0.607, 0.717, 0.442, 0.456, 0.617, 0.585, 0.691, 0.435, 0.584, 0.735, 0.663, 0.47]
tanimoto chemical similarity: [0.23, 0.357, 0.207, 0.249, 0.118, 0.282, 0.221, 0.246, 0.226, 0.212, 0.211, 0.276, 0.098, 0.259, 0.2, 0.322, 0.333, 0.212, 0.155, 0.374]
Encoded Molecule 9:
(non-aligned) shape scores: [0.487, 0.689, 0.618, 0.622, 0.66, 0.728, 0.522, 0.687, 0.667, 0.535, 0.687, 0.594, 0.802, 0.443, 0.798, 0.645, 0.555, 0.898, 0.72, 0.683]
tanimoto chemical similarity: [0.256, 0.184, 0.276, 0.281, 0.336, 0.34, 0.328, 0.303, 0.33, 0.338, 0.42, 0.337, 0.52, 0.2, 0.336, 0.331, 0.178, 0.479, 0.441, 0.326]
Encoded Molecule 10:
(non-aligned) shape scores: [0.47, 0.489, 0.557, 0.726, 0.644, 0.705, 0.43, 0.488, 0.543, 0.668, 0.402, 0.394, 0.463, 0.428, 0.388, 0.708, 0.519, 0.745, 0.265, 0.63]
tanimoto chemical similarity: [0.444, 0.347, 0.389, 0.285, 0.231, 0.255, 0.225, 0.279, 0.164, 0.325, 0.213, 0.228, 0.202, 0.21, 0.254, 0.279, 0.232, 0.352, 0.234, 0.26]

因为我们没有 OpenEye，所以只能使用 shaep 代替 OpenEye 计算生成分子与参考分子之间的形状相似性。然后提取最为形状相似的分子。

首先，需要将参考分子和生成的分子保存成sdf。保存目录为 ./test文件夹，需要先创建。

from rdkit import Chem# 保存参考分子
sdf_writer = Chem.SDWriter('./test/ref.sdf')
for i, mol in enumerate(reference_mols_list):sdf_writer.write(mol)
sdf_writer.close()# 保存生成的分子
sdf_writer = Chem.SDWriter('./test/gen.sdf')
for i, mol in enumerate(unaligned_mols_list):for k, mol_i in enumerate(mol):sdf_writer.write(mol_i)
sdf_writer.close()

使用 shaep 计算形状相似性，参考分子作为查询分子，生成分子作为分子库

结果保存在 similarity-all.txt 文件中。注意：下面的命令不能重复运行，原生成的文件需要先删除

! ../shaep -q ./test/ref.sdf  ./test/gen.sdf  ./test/similarity-all.txt --onlyshape

运行输出：

ShaEP version 1.3.1.f5e1226987eb27f7ad273010d5ca3cb8beb9fe02 64-bit build Feb 22 2020 15:04:35
Copyright (c) 2007-2020 Mikko J. Vainio and J. Santeri Puranen.
Copyright (c) 2010 Visipoint Ltd. www.visipoint.fi
All rights reserved.1 molecule and 200 structures processed in 00:00:03.130000 wall,
00:00:02.660000 user,
00:00:00 system,
00:00:02.660000 total CPU time.

使用如下函数，解析./test/similarity-all.txt，提取各参考分子最相似的数据库分子

文本解析函数：

import numpy as npdef get_similarity(path='./test/similarity-3.txt', ref_num=10, lib_num=20):with open(path, 'r') as f:lines = f.readlines()# 提取矩阵data = lines[1].strip().split('\t')[6+ref_num:]# 相似度数值化data = [float(i) if i != 'nan' else 0 for i in data ]# 矩阵data = np.array(data)data.resize(lib_num, ref_num)return data

保存每个参考分子对应的最相似分子

results = get_similarity(path='./test/similarity-all.txt', ref_num=10, lib_num=200)
index = np.argmax(results, axis=0)
index

对应的序号是：

array([  7,  34,  53,  62,  99, 116, 198, 141, 177, 197])

提取最相似的分子

max_align_mol = []# 生成分子组成库
unaligned_mols_list_ = [j for i in unaligned_mols_list for j in i]# 提取分子
for i, num in enumerate(index):sdf_writer = Chem.SDWriter('./test/best_align_{}_{}.sdf'.format(i, num))sdf_writer.write(unaligned_mols_list_[num])sdf_writer.close()max_align_mol.append(unaligned_mols_list_[num])

可视化结果，2D分子结构

from IPython.core.display import Image, display
from rdkit.Chem import Drawdef get_2D_mol(mol):mol_2D = deepcopy(mol)rdkit.Chem.rdDepictor.Compute2DCoords(mol_2D)return mol_2Dfor i in range(len(max_align_mol)):print('{}: encoded, decoded'.format(i))print('aligned shape similarity:', results[index[i], i])display(Draw.MolsToGridImage([get_2D_mol(reference_mols_list[i]), get_2D_mol(max_align_mol[i])]))

输出示例，参考分子1（左），最相似的分子右：

参考分子1为目标生成的分子示例：

以下为一些生成的3D分子示例，其中，绿色为参考分子，cayan色为AQUID模型生成的分子，紫色为经过对齐以后的分子。可以看到形状匹配的非常好。

参考分子 1， aligned shape similarity: 0.920819

参考分子2，aligned shape similarity: 0.918124

参考分子3， aligned shape similarity: 0.905439

以下是作者的源代码，使用openeye计算生成分子与原分子的形状相似性。（缺少linsece 无法运行）

计算相似性：

# Note that we can safely ignore these aromaticity warnings.from utils.openeye_utils import *all_shape_similarity = []
for r, ref in enumerate(reference_mols_list):ROCS_out = ROCS_shape_overlap(unaligned_mols_list[r], ref)aligned_shape_similarity = [ROCS_out[i][1] for i in range(len(ROCS_out))]all_shape_similarity.append(aligned_shape_similarity)

选择形状相似性最好的分子：

best_generated_indices = [np.argmax(s) for s in all_shape_similarity]
best_generated_mols = [unaligned_mols_list[i][best_generated_indices[i]] for i in range(len(best_generated_indices))]

可视化，展示最相似的分子：

from IPython.core.display import Image, display
from rdkit.Chem import Drawdef get_2D_mol(mol):mol_2D = deepcopy(mol)rdkit.Chem.rdDepictor.Compute2DCoords(mol_2D)return mol_2D# Displaying most shape-similar molecules
for i in range(len(best_generated_mols)):print('encoded, decoded')print('aligned shape similarity:', all_shape_similarity[i][best_generated_indices[i]])display(Draw.MolsToGridImage([get_2D_mol(reference_mols_list[i]), get_2D_mol(best_generated_mols[i])]))