这个教程我们将教你使用机器学习模型和分子拼接方法来预测蛋白—配体复合物的结合能。记得配体是一些与蛋白结合（通常非共价键合）的小分子。分子拼接进行几何计算以找到小分子与蛋白结合口袋（即，蛋白质的局部有个沟，小分子可以停在那里）相互作用的结合位点。蛋白质的结构可以通过Cryo-EM或X-ray晶体学技术实验确定。这对基于结构的药物发现来说是非常强大的工具。对于分子拼接的更多信息，可以参阅 AutoDock Vina paper 和 deepchem.dock文档。有许多图形用户接口和命令行接口(如AutoDock)来进行分子拼接。这里我们展示用DeepChem编程的方法来进行拼接，它可以自动化并易于集成到机器学习管线中。

随着你学习本教程，你会做如下事情：

1. 加载蛋白配体复合物数据集(PDBbind)。
2. 进行分子拼接编程。
3. 用相互作用指纹特征化蛋白-配体复合物。
4. 拟合随机森林模型并预测结合亲和力。

为了开始本教程，我们使用简单的gzipped格式预处理的数据集文件。每行是一个分子系纺，每列是该系统的不同的信息。例如，本例中，每行反映蛋白配体复合物，有如下的列：唯一的复合物标识，配体的SMILES字串，复合物中配体与蛋白的结合亲和力（Ki），所有蛋白的PDB文件列表，所有配体的文件列表。

蛋白质—配体复合物数据

进行拼接时可视化蛋白和配体是很有用的。不幸的是，Google Colab现在还不支持我们用来可视化的Jupyter widgets。需要在你的本地机器安装 MDTraj 和 nglview来可视化我们工作的蛋白配体复合物。

In [ ]:

!pip install -q mdtraj nglview

# !jupyter-nbextension enable nglview --py --sys-prefix  # for jupyter notebook

# !jupyter labextension install  nglview-js-widgets  # for jupyter lab

In [ ]:

import os

import numpy as np

import pandas as pd

import tempfile

from rdkit import Chem

from rdkit.Chem import AllChem

import deepchem as dc

from deepchem.utils import download_url, load_from_disk

为子展示拼接过程，我们使用含有配体的SMILES字串的CSV文件以及来自PDBbind的配体和蛋白靶点的PDB文件。我们也展示如何下载和特征化PDBbind来训练模型。

In [ ]:

data_dir = dc.utils.get_data_dir()

dataset_file = os.path.join(data_dir, "pdbbind_core_df.csv.gz")

if not os.path.exists(dataset_file):

    print('File does not exist. Downloading file...')

    download_url("https://s3-us-west-1.amazonaws.com/deepchem.io/datasets/pdbbind_core_df.csv.gz")

    print('File downloaded...')

raw_dataset = load_from_disk(dataset_file)

raw_dataset = raw_dataset[['pdb_id', 'smiles', 'label']]

File does not exist. Downloading file...

File downloaded...

Let's see what raw_dataset looks like:

In [ ]:

raw_dataset.head(2)

Out[ ]:

	pdb-id	smiles	label
0	2d3u	CC1CCCCC1S(O)(O)NC1CC(C2CCC(CN)CC2)SC1C(O)O	692
1	3cyx	CC(C)(C)NC(O)C1CC2CCCCC2C[NH+]1CC(O)C(CC1CCCCC...	800

修改PDB文件

下一步，我们获得一些PDB蛋白质文件来可视化和拼接。我们使用 raw_dataset的PDB IDs并从用 pdbfixer直接从Protein Data Bank 下载pdb文件。我们也用 RDKit清洗一下结构。这可以保证蛋白和配体文件的问题正确（非标准化残差，化学验证，等）。可以自由的修改单元格和pdbids来考虑新的蛋白—配体复合物。这里我们要注意PDB文件是复杂的并需要人为判断来准备进行拼接的蛋白结构。DeepChem包含大量的拼接工具来帮助你准备蛋白质文件，但是应该在拼接前检结果。

In [ ]:

from simtk.openmm.app import PDBFile

from pdbfixer import PDBFixer

from deepchem.utils.vina_utils import prepare_inputs

In [ ]:

# consider one protein-ligand complex for visualization

pdbid = raw_dataset['pdb_id'].iloc[1]

ligand = raw_dataset['smiles'].iloc[1]

In [ ]:

%%time

fixer = PDBFixer(pdbid=pdbid)

PDBFile.writeFile(fixer.topology, fixer.positions, open('%s.pdb' % (pdbid), 'w'))

p, m = None, None

# fix protein, optimize ligand geometry, and sanitize molecules

try:

    p, m = prepare_inputs('%s.pdb' % (pdbid), ligand)

except:

    print('%s failed PDB fixing' % (pdbid)) 

if p and m:  # protein and molecule are readable by RDKit

    print(pdbid, p.GetNumAtoms())

    Chem.rdmolfiles.MolToPDBFile(p, '%s.pdb' % (pdbid))

    Chem.rdmolfiles.MolToPDBFile(m, 'ligand_%s.pdb' % (pdbid))

3cyx 1415

CPU times: user 2.43 s, sys: 249 ms, total: 2.68 s

Wall time: 5.98 s

可视化

如果你没使用Colab，你可以扩展这些单元格并用MDTraj 和nglview来可视化蛋白质和配体。

In [ ]:

import mdtraj as md

import nglview

from IPython.display import display, Image

Let's take a look at the first protein ligand pair in our dataset:

In [ ]:

protein_mdtraj = md.load_pdb('3cyx.pdb')

ligand_mdtraj = md.load_pdb('ligand_3cyx.pdb')

我们使用便利的函数 nglview.show_mdtraj来可视化蛋白和配体。注意这你只有去掉上面的单元格的注释符才能工作，安装nglview，允许必要的notebook扩展。

In [ ]:

v = nglview.show_mdtraj(ligand_mdtraj)

In [ ]:

display(v)  # interactive view outside Colab

NGLWidget()

现在我们了解配体看起来像什么，来看一下我们的蛋白。

In [ ]:

view = nglview.show_mdtraj(protein_mdtraj)

display(view)  # interactive view outside Colab

NGLWidget()

分子拼接

现在我们有了数据，可视化工具，并运行了，我们看一下是否可以使用分子拼接来评估配体和蛋白的结合亲和力。

设置拼接工作分三步，你要用不同的设置进行实验。我们要指明的三件事情是：1）如何识别靶蛋白中的结合口袋；2）如何产生配体在结合口袋中的位点（几何构像）；3）如何给位置评分。记住我们的目的是要识别能与靶蛋白强烈结合的候选配体，通过分值来反映。

DeepChem有内建的方法来识别蛋白质的结合口袋。它基于convex hull方法。该方法由产生蛋白结构周围的3D polyhedron (convex hull)并识别最近convex hull的蛋白质的表面原子。考虑一些生物化学特性，因此该方法不是纯几何的。这的优点是计算成本低并很适合我们的目的。

In [ ]:

finder = dc.dock.binding_pocket.ConvexHullPocketFinder()

pockets = finder.find_pockets('3cyx.pdb')

len(pockets)  # number of identified pockets

Out[ ]:

37

姿势的产生非常复杂。幸运的是，使用DeepChem的姿势生成器将安装AutoDock Vina引擎，让我们可以快速度的产生姿势。

In [ ]:

vpg = dc.dock.pose_generation.VinaPoseGenerator()

    我们可以从deepchem.dock.pose_scoring指明姿势分值函数，它包括脉冲和疏水反应以及氢键。Vina会关注这些，因此我们允许Vina为这一目的计算分值。

In [ ]:

!mkdir -p vina_test

In [ ]:

%%time

complexes, scores = vpg.generate_poses(molecular_complex=('3cyx.pdb', 'ligand_3cyx.pdb'),  # protein-ligand files for docking,

                                       out_dir='vina_test',

                                       generate_scores=True

                                      )

当产生姿势时， 对于num_modes我们使用默认值，所以Vina会产生9个它找到的最低能量姿势，单位为kcal/mol。

In [ ]:

scores

Out[ ]:

[-8.9, -8.8, -8.8, -8.8, -8.7, -8.6, -8.6, -8.5, -8.4]

们可以看到蛋白和配体的复合物吗？可以，但是我们要组合这些分子到一个RDkit分子。

In [ ]:

complex_mol = Chem.CombineMols(complexes[0][0], complexes[0][1])

让我们来可视化我们的复合物。我们能看到配体插入到蛋白质的口袋中。

In [ ]:

v = nglview.show_rdkit(complex_mol)

display(v)

NGLWidget()

    现在我们理解了这些过程的每一步，我们可以使用DeepChem的 Docker类组合所有这些。Docker产生一个生成器，生成器能产生定位的复合物和拼接值的元组。

In [ ]:

docker = dc.dock.docking.Docker(pose_generator=vpg)

posed_complex, score = next(docker.dock(molecular_complex=('3cyx.pdb', 'ligand_3cyx.pdb'),

                                          use_pose_generator_scores=True))

结合亲和力建模

拼接是有用的，尽管对于预测蛋白-配体结合亲和力是粗糙的工具。但是它需要一些时间，特别是对于大规模的虚拟筛选实验我们需要考虑不同的靶点和成千上成的蛋白配体。我们很自然会问，我们能训练机器学习模型来预测拼接分值吗？让我们来试一试吧。

我们将展示如何下载PDBbind数据集。我们可以使用MoleculeNet的加载器从PDBbind “提炼”集或“通用”集得到4852个蛋白配体复合物。为简单起见，我们使用约100个已处理的复合物来训练模型

接下来我们要转换我们的蛋白-配体复体物到可以被机器学习算法使用的表示形式。理想地，我们有神经蛋白-配体复合物指纹，但是DeepChem还没有经过学习的这类指纹。但是我们有调试好的手工特征化器能帮助我们完成我们的挑战。本教程我们将使用二个指纹，CircularFingerprint 和 ContactCircularFingerprint。DeepChem也有voxelizers 和grid descriptors转换包含原子排列的3D立体到指纹。这些特征化器对于理解蛋白配体复合物是有用的，因为它们允许我们翻译复合物到可以传递给机器学习算法的矢量。首先，我们产生circular指纹。这些转换小分子到片断向量。

In [ ]:

pdbids = raw_dataset['pdb_id'].values

ligand_smiles = raw_dataset['smiles'].values

In [ ]:

%%time

for (pdbid, ligand) in zip(pdbids, ligand_smiles):

  fixer = PDBFixer(url='https://files.rcsb.org/download/%s.pdb' % (pdbid))

  PDBFile.writeFile(fixer.topology, fixer.positions, open('%s.pdb' % (pdbid), 'w'))

  p, m = None, None

  # skip pdb fixing for speed

  try:

    p, m = prepare_inputs('%s.pdb' % (pdbid), ligand, replace_nonstandard_residues=False,

                          remove_heterogens=False, remove_water=False,

                          add_hydrogens=False)

  except:

    print('%s failed sanitization' % (pdbid))

  if p and m:  # protein and molecule are readable by RDKit

    Chem.rdmolfiles.MolToPDBFile(p, '%s.pdb' % (pdbid))

    Chem.rdmolfiles.MolToPDBFile(m, 'ligand_%s.pdb' % (pdbid))

1hfs failed sanitization

CPU times: user 2min 54s, sys: 1.44 s, total: 2min 56s

Wall time: 3min 17s

In [ ]:

proteins = [f for f in os.listdir('.') if len(f) == 8 and f.endswith('.pdb')]

ligands = [f for f in os.listdir('.') if f.startswith('ligand') and f.endswith('.pdb')]

我们做一些清洗来保证我们对于每个有效的蛋白有有效的配体文件。这些代码将比较配体和蛋白文件的PDB IDs并去除没有相应配体的蛋白。

In [ ]:

# Handle failed sanitizations

failures = set([f[:-4] for f in proteins]) - set([f[7:-4] for f in ligands])

for pdbid in failures:

  proteins.remove(pdbid + '.pdb')

In [ ]:

len(proteins), len(ligands)

Out[ ]:

(190, 190)

In [ ]:

pdbids = [f[:-4] for f in proteins]

small_dataset = raw_dataset[raw_dataset['pdb_id'].isin(pdbids)]

labels = small_dataset.label

In [ ]:

fp_featurizer = dc.feat.CircularFingerprint(size=2048)

In [ ]:

features = fp_featurizer.featurize([Chem.MolFromPDBFile(l) for l in ligands])

In [ ]:

dataset = dc.data.NumpyDataset(X=features, y=labels, ids=pdbids)

train_dataset, test_dataset = dc.splits.RandomSplitter().train_test_split(dataset, seed=42)

 便得的加载器dc.molnet.load_pdbbind将关注于下载和特征化pdbbind数据集。这需要很长的时间并计算，因此做这些工作的代码被注释。如果你想特征化所有的PDBbind提炼集，你可以去掉注释，喝杯咖啡。不然，你就继续我们前面建构的小数据集吧。

In [ ]:

# # Uncomment to featurize all of PDBBind's "refined" set

# pdbbind_tasks, (train_dataset, valid_dataset, test_dataset), transformers = dc.molnet.load_pdbbind(

#     featurizer=fp_featurizer, set_name="refined", reload=True,

#     data_dir='pdbbind_data', save_dir='pdbbind_data')

现在我们准备学习吧。

为了拟合deepchem模型，我们首先要实例化提供的（或用户写的）模型类。这时，我们创建了下个便利的类来打包Sci-Kit Learn里可用的任何机器学习模型以便以和deepchem互操作。为了SklearnModel，你需要(a) 任务类型, (b)模型参数，另一个下面要说明的dict类，以及(c)一个 定义你要拟合的模型类型的model_instance，这里我们用RandomForestRegressor。

In [ ]:

from sklearn.ensemble import RandomForestRegressor

from deepchem.utils.evaluate import Evaluator

import pandas as pd

In [ ]:

seed = 42 # Set a random seed to get stable results

sklearn_model = RandomForestRegressor(n_estimators=100, max_features='sqrt')

sklearn_model.random_state = seed

model = dc.models.SklearnModel(sklearn_model)

model.fit(train_dataset)

注意测试集的值提示模型没有产一有意义的输出。结果证明预测结合亲和力是困难的。本教程不是要展示如何创建最先进的模型来预测结合亲和力，而是给你工具来产生你自已的分子拼接、复合物特征化数据集，并训练模型。

In [ ]:

# use Pearson correlation so metrics are > 0

metric = dc.metrics.Metric(dc.metrics.pearson_r2_score)

evaluator = Evaluator(model, train_dataset, [])

train_r2score = evaluator.compute_model_performance([metric])

print("RF Train set R^2 %f" % (train_r2score["pearson_r2_score"]))

evaluator = Evaluator(model, test_dataset, [])

test_r2score = evaluator.compute_model_performance([metric])

print("RF Test set R^2 %f" % (test_r2score["pearson_r2_score"]))

RF Train set R^2 0.921762

RF Test set R^2 0.051527

我们使用非常小的数据集并只是简单的表示，所以测试集的性能很差也不奇怪了。

In [ ]:

# Compare predicted and true values

list(zip(model.predict(train_dataset), train_dataset.y))[:5]

Out[ ]:

[(7.445426666666658, 7.4),

 (6.414323333333337, 6.85),

 (4.788719999999994, 3.4),

 (6.795444444444452, 6.72),

 (8.88874999999999, 11.06)]

In [ ]:

list(zip(model.predict(test_dataset), test_dataset.y))[:5]

Out[ ]:

[(5.917333333333333, 4.21),

 (6.076334761904759, 8.7),

 (6.695960000000001, 6.39),

 (5.6465999999999985, 4.94),

 (5.621683333333333, 9.21)]

蛋白配体复合物的观察

上一节，我们只特征化了配体。这次，我们来看一下能否用使用结构信息的蛋白配体指纺做一些直观的事情。发了开始，我们需要再特征化数据集但是这次使用contact指纹。

In [ ]:

fp_featurizer = dc.feat.ContactCircularFingerprint(size=2048)

In [ ]:

features = fp_featurizer.featurize(zip(ligands, proteins))

dataset = dc.data.NumpyDataset(X=features, y=labels, ids=pdbids)

train_dataset, test_dataset = dc.splits.RandomSplitter().train_test_split(dataset, seed=42)

现在我们用数据集训练简单的随机森林模型。

In [ ]:

seed = 42 # Set a random seed to get stable results

sklearn_model = RandomForestRegressor(n_estimators=100, max_features='sqrt')

sklearn_model.random_state = seed

model = dc.models.SklearnModel(sklearn_model)

model.fit(train_dataset)

来看一下准确度如何！

In [ ]:

metric = dc.metrics.Metric(dc.metrics.pearson_r2_score)

evaluator = Evaluator(model, train_dataset, [])

train_r2score = evaluator.compute_model_performance([metric])

print("RF Train set R^2 %f" % (train_r2score["pearson_r2_score"]))

evaluator = Evaluator(model, test_dataset, [])

test_r2score = evaluator.compute_model_performance([metric])

print("RF Test set R^2 %f" % (test_r2score["pearson_r2_score"]))

RF Train set R^2 0.416325

RF Test set R^2 0.051535

好，看来我的准确度比只有配体的数据集要低一些。但是，蛋白配体模型还是有用的因为比纯配体模型它可以学习不同的特征。