DeLinker是一个既可以实现分子骨架跃迁又可以实现分子连接的AI方法，发表于J. Chem. Inf. Model. 2020, 60, 4, 1983–1995。链接是：https://pubs.acs.org/doi/10.1021/acs.jcim.9b01120
这里介绍DeLinker在分子骨架跃迁中应用部分的代码，以及用法。GitHub中也直接给出了骨架跃迁端到端的使用方法，可以去查看一下。
Delinker模型有一些有趣的特点：
- 不仅仅是将分子链接在一起，而是会根据两个分子片段的3D信息，例如距离与取向， 进行连接，进行骨架跃迁或者生成linker；
- 非常直接的端到端的形式；
- 提供了很多关于计算化学中分子操作的代码，例如，构象生成，形状匹配，三维结构的RMSD等计算方法，并且还是多进程形式的。
这个作者简直不要太好。代码值得仔细学习。
接下来介绍一下samples中DeLinker_scaffold_hopping.ipynb中的应用方法。由于DeLinker的代码很完善，所以可以直接修改输入部分的sdf文件即可。
1. 安装DeLinker环境
DeLinker直接给出了相应的安装环境，直接conda install一下即可。
conda env create -f DeLinker-env.yml
conda activate DeLinker-env然后从GitHub上复制相应的代码：
git clone https://github.com/fimrie/DeLinker进入samples目录。
cd DeLinker/samples接下来介绍一下samples中DeLinker_scaffold_hopping.ipynb中的应用方法。由于DeLinker的代码很完善，所以可以直接修改输入部分的sdf文件即可。
.2. 骨架跃迁代码解析
（1）导入相关的模块
import sys
sys.path.append("../")
sys.path.append("../analysis/")
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem import Draw
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem.Draw import MolDrawing, DrawingOptions
from rdkit.Chem import MolStandardize
import numpy as np
from itertools import product
from joblib import Parallel, delayed
import re
from collections import defaultdict
from IPython.display import clear_output
IPythonConsole.ipython_useSVG = True（2）基本设置
# 任务使用处理器核数
n_cores = 4
# 是否使用GPU进行Delinker分子生成
use_gpu = False（3）加载起始分子，预处理
我们选择从基于吲唑的支架 (PDB ID： 3FI3) 中构建跃迁分子。 PDB ID 3FI2 是一个成功的例子，其中的氨基吡唑核心满足研究人员的需求，即 JNK3（UniProt：P53779）的选择性比 p38（UniProt：Q16539）提高 > 1000 倍。
我们先看看分子结构，左侧是基础分子，希望得到右侧的跃迁分子。
scaff_1_path = './3FI3_ligand.sdf' #起始分子
scaff_2_path = './3FI2_ligand.sdf' #跃迁以后希望得到的分子
scaff_1_sdf = Chem.SDMolSupplier(scaff_1_path)
scaff_1_smi = Chem.MolToSmiles(scaff_1_sdf[0])
scaff_2_sdf = Chem.SDMolSupplier(scaff_2_path)
scaff_2_smi = Chem.MolToSmiles(scaff_2_sdf[0])
img = Draw.MolsToGridImage([Chem.MolFromSmiles(scaff_1_smi), Chem.MolFromSmiles(scaff_2_smi)], molsPerRow=2, subImgSize=(300, 300))
img
现在要将分子1处理成DeLinker的输入形式。为此需要确定哪一部分是需要被替换的。
starting_point_2d = Chem.Mol(scaff_1_sdf[0]) #mol对象
_ = AllChem.Compute2DCoords(starting_point_2d) #生成分子的二维坐标
example_utils.mol_with_atom_index(starting_point_2d) #为每个原子生成序号
现在要选择两条键，对分子进行分段。这里通过原子来原则，例如，选择19-21原子之间和55-16原子之间的键作为断开位置，生成两侧的片段。
atom_pair_idx_1 = [19, 21]
atom_pair_idx_2 = [5, 16]
bonds_to_break = [starting_point_2d.GetBondBetweenAtoms(x,y).GetIdx() for x,y in [atom_pair_idx_1, atom_pair_idx_2]] #找到断键的位置
fragmented_mol = Chem.FragmentOnBonds(starting_point_2d, bonds_to_break) #分子根据断键位置切成片段
_ = AllChem.Compute2DCoords(fragmented_mol)
fragmented_mol #切成三部分了
确定linker以及两头的片段
fragmentation = Chem.MolToSmiles(fragmented_mol).split('.')
fragments = []
for fragment in fragmentation:
#*代表断键位置，所以有两个*号的是linker
if len([x for x in fragment if x =="*"]) ==2:
linker=fragment
else:
fragments.append(fragment)
fragments = '.'.join(fragments) #片段用.链接在一起
linker = re.sub('[0-9]+\*', '*', linker)
fragments = re.sub('[0-9]+\*', '*', fragments)生成的linker和fragments的smiles分别是：
准备数据输入到DeLinker
计算片段之间的距离和角度
dist, ang = frag_utils.compute_distance_and_angle(scaff_1_sdf[0], linker, fragments)
fragments, dist, ang写入txt文件中：
data_path = "./scaffold_hopping_test_data.txt"
with open(data_path, 'w') as f:
f.write("%s %s %s" % (fragments, dist, ang))生成scaffold_hopping_test_data.txt文件，其内容如下：
将上述的txt中的片段转化成，smi_mol, smi_linker等mol对象；然后按照ZINC的数据格式，变成分子图的节点和边的特征，保存在文件scaffold_hopping_test.json中。"zinc"这个参数一定要这么写，只要是正常的smiles作为片段的输入，都要填写zinc,这个是文件写死的格式控制。
raw_data = read_file(data_path)
preprocess(raw_data, "zinc", "scaffold_hopping_test", True)处理过程输出：
（4）使用DeLinker来生成新的分子
如果不使用GPU：
import os
if not use_gpu:
os.environ['CUDA_VISIBLE_DEVICES'] = '-1'设置DeLinker的参数。数据的格式还是选择zinc, 生成10个分子，linker的最大原子数为10，最小原子数为9，使用的输入jsion文件是之前的scaffold_hopping_test.json，train\valid都选择该文件。输出文件名为：DeLinker_example_generation_scaffold_hopping.smi
args = defaultdict(None)
args['--dataset'] = 'zinc'
args['--config'] = '{"generation": true, \
"batch_size": 1, \
"number_of_generation_per_valid": 10, \
"min_atoms": 9, "max_atoms": 10, \
"train_file": "molecules_scaffold_hopping_test.json", \
"valid_file": "molecules_scaffold_hopping_test.json", \
"output_name": "DeLinker_example_generation_scaffold_hopping.smi"}'
args['--freeze-graph-model'] = False
args['--restore'] = '../models/pretrained_DeLinker_model.pickle' #加载模型的权重设置模型生成分子，共生成20个分子。
# Setup model and generate molecules
model = DenseGGNNChemModel(args)
model.train()
# Free up some memory 清空内存
model = ''
（5）评估生成的新的分子
首先评估生成分子的2D和3D的性质，利用的DeLinker中的过滤/打分机制。
# 加载生成的分子
generated_smiles = frag_utils.read_triples_file("./DeLinker_example_generation_scaffold_hopping.smi")
in_mols = [smi[1] for smi in generated_smiles] #原始分子
frag_mols = [smi[0] for smi in generated_smiles] #分子片段
gen_mols = [smi[2] for smi in generated_smiles] #生成的分子
#生成干净的分子片段，通过分子片段替换的方式
du = Chem.MolFromSmiles('*')
clean_frags = [Chem.MolToSmiles(Chem.RemoveHs(AllChem.ReplaceSubstructs(Chem.MolFromSmiles(smi),du,Chem.MolFromSmiles('[H]'),True)[0])) for smi in frag_mols]
clear_output(wait=True)
print("Done")例如，将[*]c1cccc(C(=O)Nc2cc(OC)c(OC)c(OC)c2)c1.[*]Nc1ccccc1F两个分子片段，生成COc1cc(NC(=O)c2ccccc2)cc(OC)c1OC.Nc1ccccc1F。实际上就是补氢操作。
查看生成分子的有效性
查看生成的分子中，有多少包含原来的片段。使用子结构查询，需要先给片段加氢，也就是之前的操作方法。分子有效率是100%。
results = []
for in_mol, frag_mol, gen_mol, clean_frag in zip(in_mols, frag_mols, gen_mols, clean_frags):
if len(Chem.MolFromSmiles(gen_mol).GetSubstructMatch(Chem.MolFromSmiles(clean_frag)))>0:
results.append([in_mol, frag_mol, gen_mol, clean_frag])
print("Number of generated SMILES: \t%d" % len(generated_smiles))
print("Number of valid SMILES: \t%d" % len(results))
print("%% Valid: \t\t\t%.2f%%" % (len(results)/len(generated_smiles)*100))
找出生成的linker
这里还用到了多线程的方式，将找到的linker填充在results里面。
# 找到生成分子中的linkers
linkers = Parallel(n_jobs=n_cores)(delayed(frag_utils.get_linker)(Chem.MolFromSmiles(m[2]), Chem.MolFromSmiles(m[3]), m[1]) \
for m in results)
# linkers标准化
for i, linker in enumerate(linkers):
if linker == "":
continue
linker = Chem.MolFromSmiles(re.sub('[0-9]+\*', '*', linker))
Chem.rdmolops.RemoveStereochemistry(linker)
linkers[i] = MolStandardize.canonicalize_tautomer_smiles(Chem.MolToSmiles(linker))
# Update results
for i in range(len(results)):
results[i].append(linkers[i])
clear_output(wait=True)
print("Done")创建结果字典，检查生成不重复的分子数量。即分子生成中的uniqueness指标。
results_dict = {}
for res in results:
#去重
if res[0]+'.'+res[1] in results_dict: # Unique identifier - starting fragments and original molecule
results_dict[res[0]+'.'+res[1]].append(tuple(res))
else:
results_dict[res[0]+'.'+res[1]] = [tuple(res)]
print("Unique molecules: %.2f%%" % (frag_utils.unique(results_dict.values())*100))这里的结果是100%：
查看分子片段有效性，符合DeLinker中分子片段的比例：
filters_2d = frag_utils.calc_filters_2d_dataset(results, pains_smarts_loc="../analysis/wehi_pains.csv", n_cores=n_cores)
results_filt = []
for res, filt in zip(results, filters_2d):
if filt[0] and filt[1] and filt[2]:
results_filt.append(res)
clear_output(wait=True)
print("Pass all 2D filters: \t\t\t\t%.2f%%" % (len(results_filt)/len(results)*100))
print("Valid and pass all 2D filters: \t\t\t%.2f%%" % (len(results_filt)/len(generated_smiles)*100))
print("Pass synthetic accessibility (SA) filter: \t%.2f%%" % (len([f for f in filters_2d if f[0]])/len(filters_2d)*100))
print("Pass ring aromaticity filter: \t\t\t%.2f%%" % (len([f for f in filters_2d if f[1]])/len(filters_2d)*100))
print("Pass SA and ring filters: \t\t\t%.2f%%" % (len([f for f in filters_2d if f[0] and f[1]])/len(filters_2d)*100))
print("Pass PAINS filters: \t\t\t\t%.2f%%" % (len([f for f in filters_2d if f[2]])/len(filters_2d)*100))
获取生成的不重复的分子，结果是6个。
print("Number molecules passing 2D filters:\t\t%d" % len(results_filt))
results_filt_unique = example_utils.unique_mols(results_filt)
print("Number unique molecules passing 2D filters:\t%d" % len(results_filt_unique))
图片: https://uploader.shimo.im/f/eFD8bbihrf8yBWWH.png!thumbnail?accessToken=eyJhbGciOiJIUzI1NiIsImtpZCI6ImRlZmF1bHQiLCJ0eXAiOiJKV1QifQ.eyJleHAiOjE2NTg2MjU4MDQsImZpbGVHVUlEIjoiTDlrQk1NRFBYZUNuOHFLNSIsImlhdCI6MTY1ODYyNTUwNCwiaXNzIjoidXBsb2FkZXJfYWNjZXNzX3Jlc291cmNlIiwidXNlcklkIjo3NDE1MjkwM30.8mvl9XoJikxhjVhgohg8devHX8lN6tv5JDEZftYpaXc
5.2 3D分析
这里的目的是，为每一个生成的分子找到一个与之前分子最相似的构象！！！
用rdkit将上述不重复的分子生成构象，还是多进程的形式。rdkit比较慢，还是需要不少的时间。生成的分子构象保存在DeLinker_generated_mols_unique.sdf
_ = rdkit_conf_parallel.gen_confs([res[2] for res in results_filt_unique], "./DeLinker_generated_mols_unique.sdf",
smi_frags=[res[1] for res in results_filt_unique], numcores=n_cores, jpsettings=True)
clear_output(wait=True)
print("Done")加载分子构象，参考分子。
gen_sdfs = Chem.SDMolSupplier("./DeLinker_generated_mols_unique.sdf")
ref_mol = Chem.Mol(scaff_1_sdf[0])
img = Draw.MolsToGridImage([ref_mol], molsPerRow=1, subImgSize=(300, 300))
img
为每一种断键模式，按照断键模式，进行重新编号。
为每一种断键模式，按照断键模式，进行重新编号
# Get list of starting fragments and original molecules
used = set([])
ref_identifiers = [(res[1], res[0]) for res in results_filt if res[1]+'.'+res[0] not in used and (used.add(res[1]+'.'+res[0]) or True)]
# Get indices of compounds in SD file
start_stop_idxs = []
start = 0
errors = 0
curr_st_pt = ""
for count, gen_mol in enumerate(gen_sdfs):
try:
# Check if seen this ligand before
if gen_mol.GetProp("_Model") == str(0):
stop = count
if count != 0:
start_stop_idxs.append((start, stop))
start = int(stop) # deep copy
curr_st_pt = gen_mol.GetProp("_StartingPoint")
except:
errors += 1
continue
# Add last
start_stop_idxs.append((start, len(gen_sdfs)))
利用DeLinker自定义的SC_RDKit_frag_scores进行打分
# Calculate SC_RDKit fragments scores
names_frags = []
best_scores_frags = []
idx_best_poses_frags = []
names_frags_start_pts = []
with Parallel(n_jobs=n_cores, backend='multiprocessing') as parallel:
for i in range(-(-len(start_stop_idxs)//n_cores)):
jobs = []
for core in range(n_cores):
if i*n_cores+core < len(start_stop_idxs):
start, stop = start_stop_idxs[i*n_cores+core]
frag_smi = gen_sdfs[start].GetProp("_StartingPoint")
# Prepare jobs
gen_mols = [(Chem.Mol(gen_sdfs[idx]), Chem.Mol(ref_mol), str(frag_smi)) for idx in range(start, stop) if gen_sdfs[idx] is not None] # Test addition
jobs.append(gen_mols)
# Get SC_RDKit scores
set_scores = parallel((delayed(frag_utils.SC_RDKit_frag_scores)(gen_mols) for gen_mols in jobs))
for core, scores in enumerate(set_scores):
start, stop = start_stop_idxs[i*n_cores+core]
names_frags.append(gen_sdfs[start].GetProp("_Name"))
names_frags_start_pts.append(gen_sdfs[start].GetProp("_StartingPoint"))
best_scores_frags.append(max(scores))
idx_best_poses_frags.append((np.argmax(scores)+start, 0))
best_scores_frags_all = []
comp = list(zip(names_frags_start_pts, names_frags))
for res in results_filt:
try:
idx = comp.index((res[1], res[2]))
best_scores_frags_all.append(best_scores_frags[idx])
except:
continue计算SC_RDKit_frag_scores分数的比例：
# Print SC_RDKit Fragments results
print("Average SC_RDKit Fragments score: %.3f +- %.3f\n" % (np.mean(best_scores_frags_all), np.std(best_scores_frags_all)))
thresholds_SC_RDKit = [0.6, 0.7, 0.75, 0.8]
for thresh in thresholds_SC_RDKit:
print("SC_RDKit Fragments - Molecules above %.2f: %.2f%%" % (thresh, len([score for score in best_scores_frags_all if score >= thresh]) / len(best_score)
计算生成构象与参考分子之间构象的RMSD，找到重合最高的分子。
names_rmsd_frags = []
best_rmsd_frags = []
idx_best_rmsd_poses_frags = []
names_rmsd_frags_start_pts = []
with Parallel(n_jobs=n_cores, backend='multiprocessing') as parallel:
for i in range(-(-len(start_stop_idxs)//n_cores)):
jobs = []
for core in range(n_cores):
if i*n_cores+core < len(start_stop_idxs):
start, stop = start_stop_idxs[i*n_cores+core]
frag_smi = gen_sdfs[start].GetProp("_StartingPoint")
# Prepare jobs
gen_mols = [(Chem.Mol(gen_sdfs[idx]), Chem.Mol(ref_mol), str(frag_smi)) for idx in range(start, stop) if gen_sdfs[idx] is not None] # Test addition
jobs.append(gen_mols)
# Calculate RMSDs
set_scores = parallel((delayed(frag_utils.rmsd_frag_scores)(gen_mols) for gen_mols in jobs)) # Multiprocessing step
for core, scores in enumerate(set_scores):
start, stop = start_stop_idxs[i*n_cores+core]
names_rmsd_frags.append(gen_sdfs[start].GetProp("_Name"))
names_rmsd_frags_start_pts.append(gen_sdfs[start].GetProp("_StartingPoint"))
best_rmsd_frags.append(min(scores))
idx_best_rmsd_poses_frags.append((np.argmin(scores)+start, 0))
best_rmsd_frags_all = []
comp = list(zip(names_rmsd_frags_start_pts, names_rmsd_frags))
for res in results_filt:
try:
idx = comp.index((res[1], res[2]))
best_rmsd_frags_all.append(best_rmsd_frags[idx])
except:
continue查看分子构象重合度评分。
print("Average Fragments RMSD: %.3f +- %.3f\n" % (np.mean(best_rmsd_frags_all), np.std(best_rmsd_frags_all)))
thresholds_rmsd = [1.0, 0.75, 0.5]
for thresh in thresholds_rmsd:
print("RMSD Fragments - Molecules below %.2f: %.2f%%" % (thresh, len([score for score in best_rmsd_frags_all if score <= thresh]) / len(best_rmsd_frags_all)*100))
33%的分子的RMSD范围小于1.0A，感觉还是不错的。
（6）查看最好的分子，及其构象
有两种方式，一种是基于SC_RDKit Fragments打分排序，另一种是RMSD排序。
基于SC_RDKit Fragments打分排序：
基于SC_RDKit Fragments打分排序：
# Best by SC_RDKit Fragments
best_mols = sorted(list(zip(names_frags, best_scores_frags)), key=lambda x: x[1], reverse=True)[:10]
mols = [Chem.MolFromSmiles(m[0]) for m in best_mols]
frag_to_align = re.sub('\[\*\]', '', fragments.split('.')[0])
p = Chem.MolFromSmiles(frag_to_align)
AllChem.Compute2DCoords(p)
for m in mols: AllChem.GenerateDepictionMatching2DStructure(m,p)
Draw.MolsToGridImage(mols,
molsPerRow=3, legends=["%.2f" % m[1] for m in best_mols]
)
基于RMSD排序：
# Best by RMSD Fragments
best_mols = sorted(list(zip(names_rmsd_frags, best_rmsd_frags)), key=lambda x: x[1])[:10]
mols = [Chem.MolFromSmiles(m[0]) for m in best_mols]
frag_to_align = re.sub('\[\*\]', '', fragments.split('.')[0])
p = Chem.MolFromSmiles(frag_to_align)
AllChem.Compute2DCoords(p)
for m in mols: AllChem.GenerateDepictionMatching2DStructure(m,p)
Draw.MolsToGridImage(mols,
molsPerRow=5, legends=["%.2f" % m[1] for m in best_mols]
)
总结
1.文章提供了一种基于三维结构骨架跃迁的方法，从给出的案例上来看，效果还算可以，骨架跃迁以后的RMSD较小，形状也较为匹配；当然，也存在这生成的分子有些别扭的情况；
2.作者提供了完善的体系化的分子生成的评价方法，同时也有构象生成，分子RMSD计算等模块，并且是多线程的运行形式，非常友好。
结论是：值得尝试
657
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
