欢迎光临！ #

深度学习可以用来解决许多类型的问题，但基本工作流程通常是相同的。以下是你可以遵循的典型步骤。

选择你用来训练模型的数据集(训练集，如果没有合适的现有数据集，则建立一个新的数据集)。
建立模型。
使用训练集训练模型。
使用独立于训练集的验证集评估模型。
使用该模型对新数据进行预测。

有了DeepChem，每一个步骤都可以少到只有一两行Python代码。在本教程中，我们将通过一个基本示例演示解决现实科学问题的完整工作流。

我们要解决的问题是根据小分子的化学公式预测其溶解度。这是药物开发中一个非常重要的特性：如果一种拟用药物的溶解性不够，很可能进入患者血液的药物的量不够，从而产生不了治疗效果。我们需要的第一个东西是真实分子溶解度的数据集。DeepChem的核心组件之一，MoleculeNet，是多样化学分子数据集的集合。对于本教程，我们可以使用Delaney溶解度数据集。此数据集中的溶解度数据是以log（溶解度）表示，其中溶解度以摩尔/升为单位进行测量。

import deepchem as dc
tasks, datasets, transformers = dc.molnet.load_delaney(featurizer='GraphConv')
train_dataset, valid_dataset, test_dataset = datasets

我现在不会对这段代码解释太多。我们将在后面的教程中看到许多类似的例子。有两个细节我想让你注意一下。首先，注意传递给’ load_delaney() ‘函数的’ featurizer ‘参数。分子可以用多种方式表示。因此，我们告诉它我们想要使用哪种方式表示，或者用更专业的语言来说，如何“特征器（featurizer）”数据。其次，注意我们实际上得到了三个不同的数据集:训练集、验证集和测试集（分别对应 train_dataset, valid_dataset, test_dataset ）。在标准的深度学习工作流程中，每一个数据集都有不同的功能。

现在我们有了数据集，下一步是建立一个模型。我们将使用一种特殊的模型，称为“图卷积网络（graph convolutional network）”，简称为“graphconv”。

model = dc.models.GraphConvModel(n_tasks=1, mode='regression', dropout=0.2)

在这里我也不会过多地介绍上述代码。后面的教程将会提供更多关于“GraphConvModel”，以及DeepChem提供的其他类型的模型的信息。

我们现在需要使用训练集训练模型。我们只需给它一个数据集，然后告诉它要进行多少次训练(也就是说，要多少遍完整的使用训练集)。

model.fit(train_dataset, nb_epoch=100)

如果一切进展顺利，我们现在应该有一个完全训练好的模型了!但是这个模型靠谱吗?为了找出答案，我们必须使用验证集评估模型。我们通过选择一个评估指标并在模型上调用“evaluate()”来做到这一点。对于这个例子，让我们使用皮尔逊相关（the Pearson correlation），也称为 \(r^2\) ，作为我们的评估指标。我们可以在训练集和测试集上对它进行评估。

metric = dc.metrics.Metric(dc.metrics.pearson_r2_score)
print("Training set score:", model.evaluate(train_dataset, [metric], transformers))
print("Test set score:", model.evaluate(test_dataset, [metric], transformers))

注意，它在训练集上的得分高于测试集。对比于测试集，模型通常在训练集表现得更好。这被称为“过拟合”，这也是为什么使用独立的测试集评估模型是至关重要的原因。

我们的模型在测试集上仍然有相当不错的表现。作为比较，一个产生完全随机输出的模型的相关性为0，而一个做出完美预测的模型的相关性为1。我们的模型做得很好，所以现在我们可以用它来预测其他我们关心的分子。

因为这只是一个教程，我们没有任何其他的分子我们特别想要预测，让我们只使用测试集中的前十个分子进行预测。对于每一个分子，我们打印出其化学结构(使用SMILES字符串表示)，真实的log(溶解度)值，和预测的log(溶解度)值。

predictions = model.predict_on_batch(test_dataset.X[:10])
for molecule, test_solubility, prediction in zip(test_dataset.ids, test_dataset.y, predictions):
    print(molecule, test_solubility, prediction)

完。

处理数据集#

数据是机器学习的核心。本教程将介绍DeepChem用于存储和管理数据的“Dataset”类。它为高效地处理大量数据提供了简单但强大的工具。它还被设计成易于与其他流行的Python框架交互，如NumPy、Pandas、TensorFlow和PyTorch。

数据集的结构#

在上一个教程中，我们加载了关于化合物溶解度的Delaney数据集。现在，让我们再次加载它。

tasks, datasets, transformers = dc.molnet.load_delaney(featurizer='GraphConv')
train_dataset, valid_dataset, test_dataset = datasets

我们现在有三个Dataset对象:训练集、验证集和测试集。它们各自包含什么信息?我们可以先打印出其中一个的字符串表示形式。

print(test_dataset)

这里有很多信息，所以让我们从头开始。它以标签“DiskDataset”开始。Dataset是一个抽象类。它有几个子类，对应于存储数据的不同方式。

DiskDataset 是一个已经保存到硬盘上的数据集。数据以一种可以高效访问的方式储存在电脑上，即使数据的总量远远大于计算机的内存。
NumpyDataset 是一个存在于内存的数据集，它将所有数据保存在NumPy数组中。当操作可以完全放入内存的中小型数据集时，它是一个有用的工具。
ImageDataset 是一个更专门的类，它存储部分或所有在硬盘上的图像文件数据。在处理以图像作为输入或输出的模型时，它非常有用。

现在让我们讨论数据集的内容。每个数据集存储 样本（samples） 列表。非常粗略地说，一个样本是单个数据点。现在的情况下，每个样本都是一个分子。在其他数据集中，一个样本可能对应于一个实验测定数据、一个细胞系、一张图像或许多其他东西。对于每个样本，数据集存储以下信息。

特征，被称为“X”。这个用来作为样本输入到模型中。
标签，称为“y”。这个是我们希望模型输出的。在训练过程中，模型试图使每个样本的输出尽可能接近“y”。
权重，称为“w”。这个用来表示某些数据值比其他数据值更重要。在后面的教程中，我们将看到一些巧妙使用了权重的例子。
ID，是样本的唯一标识符。它可以是任何东西，只要它是唯一的。有时它只是一个整数索引，但在这个数据集中，ID是描述分子的SMILES字符串。

注意， X 、 y 和 w 的第一个维度的大小都是113。这意味着该数据集包含113个样本。

打印输出中列出的最后一条信息是 task_names。有些数据集包含对应于每个样本的多条信息。例如，如果一个样本代表一个分子，那么数据集可能会记录针对该分子的几个不同实验的结果。但这个数据集只有一个任务(task):“测量的 log(溶解度), 单位为摩尔/升)”。还要注意 y 和 w 都有形状(113,1)。这些数组的第二个维度通常与任务的数量相匹配。

从数据集访问数据#

有许多方法可以访问数据集中包含的数据。最简单的方法是直接访问 X ， y， w 和 ids 属性。每一个都以NumPy数组的形式返回相应的信息。

print(test_dataset.X)

这是一种非常简单的访问数据方法，但是在使用它时应该非常小心。这需要将所有样本的数据同时加载到内存中。这对于像这样的小型数据集来说没什么问题，但对于大型数据集，它很容易占用比计算机所拥有的更多的内存。

更好的方法是遍历数据集。这让它每次只加载一点数据，处理它，然后在加载下一个部分之前释放内存。你可以使用 itersamples() 方法一次遍历一个样本。

for X, y, w, id in test_dataset.itersamples():
print(X, y, w, id)

大多数深度学习模型可以同时处理一批多个样本。你可以使用 iterbatch() 来遍历每批次样本。

for X, y, w, ids in test_dataset.iterbatches(batch_size=50):
print(y.shape)

iterbatch() 在训练模型时还有其他有用的特性。例如， iterbatch(batch_size=100, epoch =10, deterministic=False) 将遍历整个数据集十次，每次都以不同的随机顺序。

数据集还可以使用TensorFlow和PyTorch的标准接口访问数据。如果要获取 tensorflow.data.Dataset ，请调用 make_tf_dataset() 。如果要获取 torch.utils.data.IterableDataset ，请调用 make_pytorch_dataset() 。有关更多细节，请参阅API文档。

最后一种访问数据的方法是 to_dataframe() 。这将数据复制到Pandas的 DataFrame 中。这需要一次性将所有数据存储在内存中，所以你应该只对小型数据集使用它。

test_dataset.to_dataframe()

创建数据集#

现在让我们谈谈如何创建自己的数据集。创建 NumpyDataset 非常简单:只需将包含数据的数组传递给构造函数。让我们创建一些随机数组，然后将它们包装在NumpyDataset中。

import numpy as np

X = np.random.random((10, 5))
y = np.random.random((10, 2))
dataset = dc.data.NumpyDataset(X=X, y=y)
print(dataset)

注意，我们没有指定权重或IDs。这些是可选的，就像 y 一样。 NumpyDataset 只要求 X 。因为我们没有给它们，它自动为我们构建 w 和 IDs 数组，将所有权重设置为1，并将IDs设置为整数索引。

dataset.to_dataframe()

如何创建 DiskDataset ？如果数据在NumPy数组中，可以调用 DiskDataset.from_numpy() 将其保存到硬盘中。由于这只是一个教程，我们将把它保存到一个临时目录。

import tempfile

with tempfile.TemporaryDirectory() as data_dir:
    disk_dataset = dc.data.DiskDataset.from_numpy(X=X, y=y, data_dir=data_dir)
    print(disk_dataset)

内存无法容纳的大型数据集怎么办?如果你在硬盘上有一些包含数以亿计分子数据的巨大文件呢?从它们创建 DiskDataset 的过程稍微复杂一些。幸运的是，DeepChem的 DataLoader 框架可以为你自动完成大部分工作。这是一个大的主题，所以我们将在后面的教程中讨论。

完。

对 MoleculeNet 的介绍#

[2] https://deepchem.readthedocs.io/en/latest/moleculenet.html

Jupyter Notebook 中文翻译版

DeepChem最强大的功能之一是它自带“电池”，也就是说，它自带数据集。DeepChem开发者社区维护MoleculeNet[1]数据集套件，它包含了大量不同的科学数据集，用于机器学习应用。最初的MoleculeNet套件有17个主要关注分子性质的数据集。在过去的几年里，MoleculeNet已经发展成为一个更广泛的科学数据集集合，以促进科学机器学习工具的广泛使用和发展。

这些数据集与DeepChem套件的其他部分集成在一起，因此你可以通过 dc.molnet 子模块中的函数方便地访问这些数据集。在学习本系列教程的过程中，你已经看到了这些加载器（loaders）的一些示例。MoleculeNet套件的完整文档可在我们的文档[2]中找到。

[1] Wu, Zhenqin, et al. “MoleculeNet: a benchmark for molecular machine learning.” Chemical science 9.2 (2018): 513-530.

MoleculeNet的概述#

在上两个教程中，我们加载了分子溶解度——Delaney数据集。让我们再加载一次。

import deepchem as dc

tasks, datasets, transformers = dc.molnet.load_delaney(featurizer='GraphConv')
train_dataset, valid_dataset, test_dataset = datasets

注意我们调用的加载器函数 dc.molnet.load_delaney 在 dc.molnet 的子模块。让我们看一看可供我们使用的加载器函数的完整名单。

[method for method in dir(dc.molnet) if "load_" in method ]

MoleculeNet加载器由DeepChem社区积极维护，我们致力于向集合中添加新的数据集。让我们看看今天在MoleculeNet中有多少数据集。

len([method for method in dir(dc.molnet) if "load_" in method ])

MoleculeNet数据集类别#

MoleculeNet 中有很多不同的数据集。让我们快速概述一下可用的不同类型的数据集。我们将把数据集分成不同的类别，并列出属于这些类别的加载器。更多关于这些数据集的细节可以在 https://deepchem.readthedocs.io/en/latest/moleculenet.html 上找到。最初的MoleculeNet论文[1]提供了标记为 V1 的数据集的详细信息。所有剩下的数据集都是 V2 ，在论文中没有记录。

量子力学数据集#

MoleculeNet的量子力学数据集包含各种量子力学特性预测任务。目前的量子力学数据集包括QM7、QM7b、QM8、QM9。相关的加载器是：

物理化学数据集#

物理化学数据集包含预测分子的各种物理性质的各种任务。

dc.molnet.load_delaney ：V1 这个数据集在原论文中也称为ESOL。
dc.molnet.load_sampl ：V1 这个数据集在原论文中也称为FreeSolv。
dc.molnet.load_lipo ：V1 这个数据集在原文中也被称为Lipophilicity。
dc.molnet.load_thermosol ：V2
dc.molnet.load_hppb ：V2
dc.molnet.load_hopv ：V2 该数据集来自最近出版的一篇文章 [3]

化学反应数据集#

这些数据集包含了用于计算逆向合成/正向合成的化学反应数据集。

dc.molnet.load_uspto

生物化学/生物物理数据集#

这些数据集是从各种生化/生物物理数据集中提取的，这些数据集包括比如化合物与蛋白质的结合亲和力等数据。

dc.molnet.load_pcba ：V1
dc.molnet.load_nci ：V2
dc.molnet.load_muv ：V1
dc.molnet.load_hiv ：V1
dc.molnet.load_ppb ：V2。
dc.molnet.load_bace_classification ：V1 这个加载器加载原MoleculeNet论文中的分类BACE数据集的程序。
dc.molnet.load_bace_regression ：V1 这个加载器加载原MoleculeNet论文中的回归分析BACE数据集的程序。
dc.molnet.load_kaggle ：V2 该数据集来自Merck的药物发现kaggle竞赛内容，在[4]中进行了描述。
dc.molnet.load_factors ：V2 该数据集来自[4]。
dc.molnet.load_uv ：V2 该数据集来自[4]。
dc.molnet.load_kinase ：V2 该数据集来自[4]。

分子目录数据集#

这些数据集提供的分子数据集除了原始的SMILES公式或结构外没有相关的性质。这种类型的数据集对于建立生成模型任务非常有用。

生理学数据集#

这些数据集包含关于化合物如何与人类患者相互作用的生理特性的实验数据。

结构生物学数据集#

这些数据集包含大分子的三维结构和相关的性质。

dc.molnet.load_pdbbind ：V1

显微术数据集#

这些数据集包含显微术图像数据集，通常是细胞系。这些数据集并没有出现在最初的MoleculeNet论文中。

材料属性数据集#

这些数据集包含关于各种材料的性能的数据。

[3] Lopez, Steven A., et al. “The Harvard organic photovoltaic dataset.” Scientific data 3.1 (2016): 1-7.

[4] Ramsundar, Bharath, et al. “Is multitask deep learning practical for pharma?.” Journal of chemical information and modeling 57.8 (2017): 2068-2076.

MoleculeNet加载器（Loaders）的解释#

所有的MoleculeNet加载器函数都采用 dc.molnet.load_X 的形式。加载器函数返回一个元组 (任务，数据集，转换器)(tasks, datasets, transformers) 。让我们遍历每个返回值并解释我们得到了什么:

任务（tasks）：这是一个任务名称列表。MoleculeNet中的许多数据集都是“多任务”的。也就是说，一个给定的数据点有多个与之相关的标签。这些对应于与这个数据点相关的不同测量值或值。
数据集（datasets）：这是一个元组包含三个 dc.data.Dataset 对象 (训练，验证，测试) 。这些对应于这个MoleculeNet数据集的训练、验证和测试集。
转换器（transformers）：这是一个在处理期间应用于此数据集的 dc.trans.Transformer 对象列表。

这有点抽象，所以让我们看看我们上面调用的 dc.molnet.load_delaney 函数的这些返回值。让我们从 任务（tasks） 开始。

print(tasks)

我们在这个数据集中有一个任务，它对应于测量的 log(溶解度)，单位为mol/L。现在让我们来看看 数据集（datasets） ：

print(datasets)

正如我们前面提到的，我们看到 datassets 是一个包含3个数据集的元组。我们把它们分开。

train, valid, test = datasets
print(train)
print(valid)
print(test)

让我们来看看 train 数据集中的一个数据点。

print(train.X[0])

注意，这是一个由 dc.feat.ConvMolFeaturizer 生成的 dc.feat.mol_graphs.ConvMol 对象。稍后我们将更多地讨论如何选择特征化（featurization）。最后让我们来看看 transformers ：

print(transformers)

我们看到一个转换器(transformer)被应用了， dc.trans.NormalizationTransformer 。

在阅读完这篇描述之后，你可能想知道在底层做了哪些选择。正如我们之前简要提到的，可以使用不同的“featurizer”来处理数据集。在这儿，我们能选择如何特征化吗?此外，如何将源数据集分割为训练/验证/测试三个不同的数据集?

你可以使用 featurizer 和 splitter 关键字参数并传入不同的字符串。“featurizer”通常可能的选择是“ECFP”，“GraphConv”，“Weave”和“smiles2img”，对应于 dc.feat.CircularFingerprint 、 dc.feat.ConvMolFeaturizer 、 dc.feat.WeaveFeaturizer 和 dc.feat.SmilesToImage 。splitter的常见可能选项是“None”，“index”，“random”，“scaffold”和“stratified”，对应于no split， dc.splits.IndexSplitter , dc.splits.RandomSplitter , dc.splits.SingletaskStratifiedSplitter 。我们还没有讨论分离器，但直观地说，它们是一种基于不同标准划分数据集的方法。我们将在以后的教程中详细介绍。

除了字符串，你还可以传入任何 Featurizer 或 Splitter 对象。这是非常有用的，例如，Featurizer的构造参数可以被用来定制它的行为。

tasks, datasets, transformers = dc.molnet.load_delaney(featurizer="ECFP", splitter="scaffold")
(train, valid, test) = datasets
print(train)
print(train.X[0])

注意，与前面的调用不同，我们有了由 dc.feat.CircularFingerprint 生成的numpy数组。而不是 dc.feat.ConvMolFeaturizer 生成的 ConvMol 对象。

自己试试吧。尝试调用MoleculeNet来加载一些其他的数据集，并使用不同的 featurizer/splitter 选项进行实验，看看会发生什么!

分子指纹#

分子可以用多种方式表示。本教程介绍一种称为“分子指纹”的表示类型。这是一种非常简单的表示，通常适用于小的类药物分子。

我们现在导入 deepchem 包来玩。

import deepchem as dc

什么是指纹?#

深度学习模型几乎总是以数字数组作为输入。如果我们想用它们来处理分子，我们需要以某种方式将每个分子表示为一个或多个数字数组。

许多(但不是所有)类型的模型要求它们的输入具有固定的大小。这对分子来说是一个挑战，因为不同的分子有不同数量的原子。如果我们想使用这些类型的模型，我们需要用固定大小的数组来表示大小不一样的分子。

指纹（fingerprints）的设计就是为了解决这些问题。指纹是一个固定长度的数组，其中不同的数位表示分子中存在不同的特征。如果两个分子有相似的指纹，这表明它们包含许多相同的特征，因此很可能具有相似的化学成分。

DeepChem支持一种特殊类型的指纹，称为“扩展连通性指纹（Extended Connectivity Fingerprint）”，简称“ECFP”。它们有时也被称为“圆形指纹”。ECFP算法首先根据原子的直接性质和化学键对其进行分类。每一个独特的图案都是一个特征。例如，“碳原子连接两个氢原子和两个重原子”将是一个特征，对于包含该特征的任何分子，指纹的一个特定数位被设置为1。然后，它通过观察更大的圆形邻域来迭代地识别新特征。一个特定的特征与另外两个特定的特征结合成为一个更高层次的特征，任何包含它的分子的分子指纹中对应的数位会被设置。这将持续进行固定数量的迭代，通常是两次。

让我们来看一个使用ECFP的数据集。

tasks, datasets, transformers = dc.molnet.load_tox21(featurizer='ECFP')
train_dataset, valid_dataset, test_dataset = datasets
print(train_dataset)

特征数组 X 的形状为 (6264,1024)。这意味着在训练集中有6264个样本。每一个都由一个长度为1024的指纹表示。还要注意标签数组 y 的形状为(6264,12):这是一个多任务数据集。Tox21包含有关分子毒性的信息。12种不同的检测方法被用来寻找毒性迹象。数据集记录了所有12个试验的结果，每一个都作为不同的任务。

让我们再看一下权重（weights）数组。

train_dataset.w

注意，有些元素是0。这些权重值用于指示丢失的数据。并不是所有的试验都是对每个分子进行的。将样本或样本/任务（tasks）对的权重设置为0将导致在拟合和评估期间忽略它。它不会对损失函数或其他指标产生影响。

大多数其他权重都接近1，但不完全是1。这样做是为了平衡每个任务中正（positive）和负（negative）样本的总权重。在训练模型时，我们希望12个任务（tasks）中的每个任务都能做出同等的贡献，并且在每个任务上，我们希望对正样本和负样本给予同等的权重。否则，模型可能只知道大多数训练样本是无毒的，因此会倾向于识别其他分子是无毒的。

基于指纹训练模型#

让我们训练一个模型。在之前的教程中，我们使用了 GraphConvModel ，这是一个相当复杂的架构，需要一组复杂的输入（inputs）。因为指纹非常简单，只是一个固定长度的数组，我们可以使用一种更简单的模型类型。

model = dc.models.MultitaskClassifier(n_tasks=12, n_features=1024, layer_sizes=[1000])

MultitaskClassifier 是一个简单的全连接层堆栈。在这个例子中，我们告诉它使用一个具有1000个神经元的隐藏层。我们还告诉它，每个输入将有1024个特征，并且它应该为12个不同的任务产生预测。

为什么不为每个任务训练一个单独的模型呢？我们可以这样做，但事实证明，训练一个多任务模型通常效果更好。我们将在后面的教程中看到其中一个例子。

让我们对模型进行训练和评估。

import numpy as np

model.fit(train_dataset, nb_epoch=10)
metric = dc.metrics.Metric(dc.metrics.roc_auc_score)
print('training set score:', model.evaluate(train_dataset, [metric], transformers))
print('test set score:', model.evaluate(test_dataset, [metric], transformers))

对于如此简单的模型和特征化来说，这是不错的结果。更复杂的模型在这个数据集上的表现稍好一些，但并不是非常好。

使用 TensorFlow 和 PyTorch 创建模型#

在之前的教程中，我们使用的是 DeepChem 提供的标准模型。这对于许多应用来说都没问题，但是迟早你会希望使用你自己定义的体系结构创建一个全新的模型。DeepChem提供了与 TensorFlow (Keras) 和 PyTorch 的集成，所以你可以在这两个框架的模型中使用它。

实际上，在DeepChem中使用TensorFlow或PyTorch模型时，你可以采用两种不同的方法。这取决于你是想使用TensorFlow/PyTorch APIs 还是DeepChem APIs 来训练和评估你的模型。对于前一种情况，DeepChem的 Dataset 类有一些方法可以方便地将其与其他框架一起使用。 make_tf_dataset() 返回一个遍历数据的 tensorflow.data.Dataset 对象。 make_pytorch_dataset() 返回一个遍历数据的 torch.utils.data.IterableDataset 对象。这让你可以使用 DeepChem 的数据集（datasets）、加载器（loaders）、特征器（featurizers）、变换器（transformers）、拆分器（splitters）等，并轻松将它们集成到你现有的 TensorFlow 或 PyTorch 代码中。

但 DeepChem 还提供了许多其他有用的功能。另一种让你使用这些功能的方法是将你的模型包装在一个 DeepChem 的 Model 对象中。让我们看看如何做到这一点。

KerasModel#

KerasModel 是DeepChem Model 类的子类。它充当 tensorflow.keras.Model 的包装器。让我们看一个使用它的例子。对于本例，我们创建了一个由两个密集层组成的简单顺序（sequential）模型。

import deepchem as dc
import tensorflow as tf

keras_model = tf.keras.Sequential([
    tf.keras.layers.Dense(1000, activation='relu'),
    tf.keras.layers.Dropout(rate=0.5),
    tf.keras.layers.Dense(1)
])
model = dc.models.KerasModel(keras_model, dc.models.losses.L2Loss())

对于本例，我们使用了Keras的 Sequential 类。我们的模型由一个具有 ReLU 激活的密集层、50% 的 dropout 提供正则化和最后一个产生标量输出的层组成。我们还需要指定在训练模型时使用的损失函数，在本例中 \(L^2\) 损失。我们现在可以训练和评估该模型，就像我们对任何其他 DeepChem 模型一样。例如，让我们加载 Delaney 溶解度数据集。我们的模型如何基于分子的扩展连通性指纹 (ECFPs) 来预测分子的溶解性?

tasks, datasets, transformers = dc.molnet.load_delaney(featurizer='ECFP', splitter='random')
train_dataset, valid_dataset, test_dataset = datasets
model.fit(train_dataset, nb_epoch=50)
metric = dc.metrics.Metric(dc.metrics.pearson_r2_score)
print('training set score:', model.evaluate(train_dataset, [metric]))
print('test set score:', model.evaluate(test_dataset, [metric]))

TorchModel#

TorchModel 的工作原理与 KerasModel 类似，只不过它包装了一个 torch.nn.Module。让我们使用PyTorch来创建另一个模型，就像前面的模型一样，并用相同的数据训练它。

import torch

pytorch_model = torch.nn.Sequential(
    torch.nn.Linear(1024, 1000),
    torch.nn.ReLU(),
    torch.nn.Dropout(0.5),
    torch.nn.Linear(1000, 1)
)
model = dc.models.TorchModel(pytorch_model, dc.models.losses.L2Loss())

model.fit(train_dataset, nb_epoch=50)
print('training set score:', model.evaluate(train_dataset, [metric]))
print('test set score:', model.evaluate(test_dataset, [metric]))

损失的计算#

现在让我们看一个更高级的例子。在上述模型中，损失是直接从模型的输出计算出来的。这通常是可以的，但并非总是如此。考虑一个输出概率分布的分类模型。虽然从概率中计算损失是可能的，但从 logits 中计算损失在数值上更稳定。

为此，我们创建一个返回多个输出的模型，包括概率和 logits。 KerasModel 和 TorchModel 让你指定一个“输出类型（output_types）”列表。如果一个特定的输出具有 ‘prediction’ 类型，这意味着它是一个正常的输出，在调用 predict() 时应该返回。如果它有 ‘loss’ 类型，这意味着它应该传递给损失函数，而不是正常的输出。

顺序模型不允许多个输出，因此我们使用子类化样式模型（subclassing style model）。

class ClassificationModel(tf.keras.Model):

    def __init__(self):
        super(ClassificationModel, self).__init__()
        self.dense1 = tf.keras.layers.Dense(1000, activation='relu')
        self.dense2 = tf.keras.layers.Dense(1)

    def call(self, inputs, training=False):
        y = self.dense1(inputs)
        if training:
            y = tf.nn.dropout(y, 0.5)
        logits = self.dense2(y)
        output = tf.nn.sigmoid(logits)
        return output, logits

keras_model = ClassificationModel()
output_types = ['prediction', 'loss']
model = dc.models.KerasModel(keras_model, dc.models.losses.SigmoidCrossEntropy(), output_types=output_types)

我们可以在 BACE 数据集中训练我们的模型。这是一个二元分类任务，试图预测一个分子是否会抑制BACE-1酶。

tasks, datasets, transformers = dc.molnet.load_bace_classification(feturizer='ECFP', splitter='scaffold')
train_dataset, valid_dataset, test_dataset = datasets
model.fit(train_dataset, nb_epoch=100)
metric = dc.metrics.Metric(dc.metrics.roc_auc_score)
print('training set score:', model.evaluate(train_dataset, [metric]))
print('test set score:', model.evaluate(test_dataset, [metric]))

类似地，我们将创建一个自定义的分类器模型类来与 TorchModel 一起使用。理由跟与上面的 KerasModel 相似，自定义模型允许轻松得到第二个密集层的未缩放输出(Tensorflow 中的 logits)。自定义类允许定义如何向前传递；在最终的sigmoid被应用产生预测之前得到 logits。

最后，用一个需要概率和 logits 的 ClassificationModel 的实例与一个损失函数搭配生成一个 TorchModel 的实例进行训练。

class ClassificationModel(torch.nn.Module):

    def __init__(self):
        super(ClassificationModel, self).__init__()
        self.dense1 = torch.nn.Linear(1024, 1000)
        self.dense2 = torch.nn.Linear(1000, 1)

    def forward(self, inputs):
        y = torch.nn.functional.relu( self.dense1(inputs) )
        y = torch.nn.functional.dropout(y, p=0.5, training=self.training)
        logits = self.dense2(y)
        output = torch.sigmoid(logits)
        return output, logits

torch_model = ClassificationModel()
output_types = ['prediction', 'loss']
model = dc.models.TorchModel(torch_model, dc.models.losses.SigmoidCrossEntropy(), output_types=output_types)

我们将使用相同的 BACE 数据集。和以前一样，该模型将尝试进行二元分类任务，试图预测一个分子是否会抑制BACE-1酶。

tasks, datasets, transformers = dc.molnet.load_bace_classification(feturizer='ECFP', splitter='scaffold')
train_dataset, valid_dataset, test_dataset = datasets
model.fit(train_dataset, nb_epoch=100)
metric = dc.metrics.Metric(dc.metrics.roc_auc_score)
print('training set score:', model.evaluate(train_dataset, [metric]))
print('test set score:', model.evaluate(test_dataset, [metric]))

其他功能#

KerasModel 和 TorchModel 有很多其他的功能。下面列一些比较重要的。

训练过程中自动保存检查点（checkpoints）。
将进度记录到控制台（console）或者传送到 TensorBoard 或 Weights & Biases 。
以 f(输出，标签，权重) 的形式定义损失函数。
使用 ValidationCallback 类来提前停止。
从已训练的模型加载参数。
估计模型输出的不确定性。
通过显著性映射（saliency mapping）识别重要特征。

通过将你自己的模型包装在 KerasModel 或 TorchModel 中，你就可以使用所有这些功能。有关它们的详细信息，请参阅API文档。

图卷积（Graph convolutions）的介绍#

在本教程中，我们将学习更多关于“图卷积”的知识。这是处理分子数据最强大的深度学习工具之一。这样做的原因是分子可以很自然地被看作图形。

https://github.com/deepchem/deepchem/blob/master/examples/tutorials/assets/basic_graphs.gif?raw=1

请注意，我们从高中开始习惯的那种标准化学图表是如何自然地将分子可视化为图表的。在本教程的其余部分中，我们将更详细地挖掘这种关系。这将让我们更深入地了解这些系统是如何工作的。

什么是图卷积？#

考虑一种通常用于处理图像的标准卷积神经网络(CNN)。输入是一个像素网格。每个像素都有一个向量，例如包含红色、绿色和蓝色通道值。数据经过一系列的卷积层。每一层都将来自像素及其邻居的数据结合起来，为像素生成一个新的数据向量。早期的层发现小规模的局部图案，而后期的层发现更大、更抽象的图案。卷积层通常与池化层交替，池化层在局部区域上执行一些操作，如筛选出最大值或最小值。

图卷积与之相似，但它们作用于图表上。它们从图表每个节点的数据向量开始(例如，该节点所代表的原子的化学性质)。卷积层和池化层将来自相互连接的节点(例如，相互连接的原子)的信息结合起来，生成一个新的数据向量作为新的节点。

训练一个GraphConvModel#

让我们使用 MoleculeNet 套件来加载 Tox21 数据集。为了以图卷积网络可以使用的方式特征化数据，我们将特征器选项设置为 ‘GraphConv’ 。MoleculeNet 调用返回一个训练集、一个验证集和一个测试集供我们使用。它还返回 tasks (任务名称的列表)和 transformers (用于预处理数据集的数据变换器的列表)。(大多数深度网络相当挑剔，需要一组数据变换器来确保训练的稳定进行。)

import deepchem as dc

tasks, datasets, transformers = dc.molnet.load_tox21(featurizer='GraphConv')
train_dataset, valid_dataset, test_dataset = datasets

现在让我们在这个数据集上训练一个图卷积网络。DeepChem 有一个名为 GraphConvModel 的类，它在底层封装了一个标准的图卷积架构，方便用户使用。让我们实例化该类的一个对象，并在数据集上训练它。

n_tasks = len(tasks)
model = dc.models.GraphConvModel(n_tasks, mode='classification')
model.fit(train_dataset, nb_epoch=50)

让我们试着评估我们训练过的模型的性能。为此，我们需要定义一个衡量标准，衡量一个模型性能的标准。 dc.metrics 是一些衡量标准的集合。对于这个数据集，使用 ROC-AUC 评分是标准的，即受试者工作特征曲线(测量精度和召回率之间的权衡)下的面积。幸运的是，DeepChem 已经提供了 ROC-AUC 评分。

为了在此衡量标准下衡量模型的性能，我们可以使用函数 model.evaluate() 。

metric = dc.metrics.Metric(dc.metrics.roc_auc_score)
print('Training set score:', model.evaluate(train_dataset, [metric], transformers))
print('Test set score:', model.evaluate(test_dataset, [metric], transformers))

结果非常好， GraphConvModel 非常容易使用。但这模型下到底发生了什么？我们可以自己构建 GraphConvModel 吗？当然可以！DeepChem为图卷积中涉及的所有计算提供Keras层。我们将使用来自 DeepChem 的以下层。

GraphConv 层：该层实现了图卷积。图卷积以一种非线性的方式将每个节点的特征向量与相邻节点的特征向量结合起来。它“混合”在一个图的局部区域中的信息。
GraphPool 层：该层对某区域内原子的特征向量进行最大池化。你可以将此层看作类似于二维卷积的最大池化层，但它对图表进行操作。
GraphGather ：许多图卷积网络对每个图标节点的特征向量进行操作。例如，对于一个分子，每个节点可能代表一个原子，网络操作代表原子的局部化学性质的特征向量。然而，最终，我们可能需要使用分子级别的特征表示。该层通过组合所有节点的特征向量来创建一个图标级别的特征向量。

除此之外，我们还将应用标准的神经网络层，如 Dense 、 BatchNormalization 和 Softmax 层。

from deepchem.models.layers import GraphConv, GraphPool, GraphGather
import tensorflow as tf
import tensorflow.keras.layers as layers

batch_size = 100

class MyGraphConvModel(tf.keras.Model):

def __init__(self):
    super(MyGraphConvModel, self).__init__()
    self.gc1 = GraphConv(128, activation_fn=tf.nn.tanh)
    self.batch_norm1 = layers.BatchNormalization()
    self.gp1 = GraphPool()

    self.gc2 = GraphConv(128, activation_fn=tf.nn.tanh)
    self.batch_norm2 = layers.BatchNormalization()
    self.gp2 = GraphPool()

    self.dense1 = layers.Dense(256, activation=tf.nn.tanh)
    self.batch_norm3 = layers.BatchNormalization()
    self.readout = GraphGather(batch_size=batch_size, activation_fn=tf.nn.tanh)

    self.dense2 = layers.Dense(n_tasks*2)
    self.logits = layers.Reshape((n_tasks, 2))
    self.softmax = layers.Softmax()

def call(self, inputs):
    gc1_output = self.gc1(inputs)
    batch_norm1_output = self.batch_norm1(gc1_output)
    gp1_output = self.gp1([batch_norm1_output] + inputs[1:])

    gc2_output = self.gc2([gp1_output] + inputs[1:])
    batch_norm2_output = self.batch_norm1(gc2_output)
    gp2_output = self.gp2([batch_norm2_output] + inputs[1:])

    dense1_output = self.dense1(gp2_output)
    batch_norm3_output = self.batch_norm3(dense1_output)
    readout_output = self.readout([batch_norm3_output] + inputs[1:])

    logits_output = self.logits(self.dense2(readout_output))
    return self.softmax(logits_output)

我们现在可以更清楚地看到正在发生的事情。有两个卷积块，每个块由一个 GraphConv 组成，然后是批正则化，然后是一个 GraphPool 来进行最大池化。最后，我们使用一个密集层、另一个批正则化、一个 GraphGather 来组合来自所有不同节点的数据，以及一个最终的密集层来生成最终输出。

现在让我们创建 DeepChem模型，它将是我们刚刚创建的 Keras 模型的包装。我们还将指定损失函数，以便模型知道最小化的目标。

model = dc.models.KerasModel(MyGraphConvModel(), loss=dc.models.losses.CategoricalCrossEntropy())

这个模型的输入是什么？图卷积需要对每个分子的完整描述，包括节点(原子)的列表，以及哪些原子彼此相连的描述。事实上，如果我们检查数据集，我们会看到特征数组包含 ConvMol 类型的Python对象。

print(test_dataset.X[0])

模型期望数字数组作为输入，而不是Python对象。我们必须将 ConvMol 对象转换为 GraphConv 、 GraphPool 和 GraphGather 层所期望的特定数组集。幸运的是， ConvMol 类包含执行此操作的代码，以及将所有分子按批次合并以创建单个数组集的代码。

下面的代码创建一个Python生成器，给定一组数据，它生成值为 Numpy 数组的输入、标签和权重列表。 atom_features 是长度为75的每个原子的特征向量。TensorFlow 要求其他输入支持小批处理。 degree_slice 是一个索引便利工具，可以方便地根据给定度数从所有分子中定位原子。 membership 表示分子跟分子中原子的隶属关系(原子 i 属于分子 membership[i] )。 deg_adjs 是一个列表，包含按原子度数分组的邻接表。要了解更多细节，请查看代码。

from deepchem.metrics import to_one_hot
from deepchem.feat.mol_graphs import ConvMol
import numpy as np

def data_generator(dataset, epochs=1):
for ind, (X_b, y_b, w_b, ids_b) in enumerate(dataset.iterbatches(batch_size, epochs,
                                                                deterministic=False, pad_batches=True)):
    multiConvMol = ConvMol.agglomerate_mols(X_b)
    inputs = [multiConvMol.get_atom_features(), multiConvMol.deg_slice, np.array(multiConvMol.membership)]
    for i in range(1, len(multiConvMol.get_deg_adjacency_lists())):
    inputs.append(multiConvMol.get_deg_adjacency_lists()[i])
    labels = [to_one_hot(y_b.flatten(), 2).reshape(-1, n_tasks, 2)]
    weights = [w_b]
    yield (inputs, labels, weights)

现在，我们可以使用 fit_generator(generator) 来训练模型，它将使用我们已经定义的生成器来训练模型。

model.fit_generator(data_generator(train_dataset, epochs=50))

现在我们已经训练了我们的图卷积网络，让我们来评估它的性能。我们必须再次使用我们定义的生成器来评估模型的性能。

print('Training set score:', model.evaluate_generator(data_generator(train_dataset), [metric], transformers))
print('Test set score:', model.evaluate_generator(data_generator(test_dataset), [metric], transformers))

成功!我们构建的模型的行为与 GraphConvModel 几乎相同。如果你希望构建你自己的定制模型，你可以遵循我们在这里提供的示例来实现这一点。我们希望很快看到你方令人兴奋的工程!

深入分子特征化#

对分子数据进行机器学习的最重要步骤之一是将数据转换为适合应用在学习算法的形式。这个过程被广泛地称为“特征化”，包括将一个分子转化为某种向量或张量。有许多不同的方法可以做到这一点，特征化方法的选择通常取决于手头的问题。我们已经看到了两种这样的方法：分子指纹和用于图卷积的 ConvMol 对象。在本教程中，我们将讨论其他一些方法。

特征器（Featurizers）#

在 DeepChem 中，将分子(或任何其他类型的输入)特征化的方法由 Featurizer 对象定义。使用特征器有三种不同的方式。

在使用 MoleculeNet 加载器函数时，你只需传递要使用的特征化方法的名称。我们已经在之前的教程中看到过这样的例子，例如 featurizer=’ECFP’ 或 featurizer=’GraphConv’ 。
你也可以创建一个特征器，并直接应用到分子。例如:

import deepchem as dc

featurizer = dc.feat.CircularFingerprint()
print(featurizer(['CC', 'CCC', 'CCO']))

在使用 DataLoader 框架创建新数据集时，可以指定用于处理数据的特征器。我们将在以后的教程中看到这一点。

我们使用丙烷( \(CH_3CH_2CH_3\) ，SMILES字符串 ‘CCC’ )作为本教程的输入。许多特征化方法使用分子的构象异构体。可以使用 deepchem.utils.conformers 中的 ConformerGenerator 类生成构象异构体。

RDKitDescriptors#

RDKitDescriptors 通过使用RDKit计算描述符列表的值来描述一个分子。这些是基本的物理和化学性质：分子量，极性表面积，氢键供体和受体的数量等。这对于预测依赖于这些高级性质而不是详细的分子结构的性质是最有用的。

特征器的本质是一组允许的描述符，可以使用 RDKitDescriptors.allowedDescriptors 来访问。此特征器使用 rdkit.Chem.Descriptors.descList 中的描述符，检查它们是否在允许的描述符列表中，并计算分子的描述符值。

让我们打印出丙烷的前十个描述符的值。

rdkit_featurizer = dc.feat.RDKitDescriptors()
features = rdkit_featurizer(['CCC'])[0]
for feature, descriptor in zip(features[:10], rdkit_featurizer.descriptors):
    print(descriptor, feature)

当然，还有比这更多的描述符。

print('The number of descriptors present is: ', len(features))

WeaveFeaturizer 和 MolGraphConvFeaturizer#

我们之前研究过图卷积，它使用 ConvMolFeaturizer 将分子转换为 ConvMol 对象。图卷积是将分子表示为图标的一类大型架构的一种特殊情况。它们的工作方式相似，但在细节上有所不同。例如，它们可能把原子或连接它们的键或这两者都用数据向量表示。他们可能会使用各种技术从前一层的数据向量计算出新的数据向量，并在最后使用各种技术计算分子级别的性质。

DeepChem 支持许多不同的基于图表的模型。其中一些需要分子被以稍微不同的方式特征化。正因为如此，有另外两个特征器称为 WeaveFeaturizer 和 MolGraphConvFeaturizer 。它们各自将分子转换为特定模型使用的不同类型的Python对象。当使用任何基于于图表的模型时，只需检查文档，看看需要使用什么特征器。

CoulombMatrix#

到目前为止，我们所研究的所有模型都只考虑了分子的内在特性:组成分子的原子列表以及连接它们的键。当处理柔性分子时，你可能还想考虑分子可以呈现的不同构象。例如，当药物分子与蛋白质结合时，结合的强度取决于原子对之间的特定相互作用。为了预测结合强度，你可能需要考虑各种可能的构象，并在进行预测时使用将这些构象考虑在内的模型。

库仑矩阵是分子构象的一种常用特征。回想一下，两个电荷之间的静电库仑相互作用与 \(q_1q_2/r\) 成正比，其中 \(q_1和q_2\) 是电荷， r 是它们之间的距离。对于一个有N原子的分子，库仑矩阵是一个 N × N 矩阵，其中每个元素给出了两个原子之间静电相互作用的强度。它包含了原子上的电荷和原子间距离的信息。更多关于函数形式的信息可以在这里找到。

为了应用这个特征器，我们首先需要分子的一组构象。我们可以使用 ConformerGenerator 类来做到这一点。它取一个 RDKit 分子，生成一组能量最小化的构象，并对其进行修剪，使其只包含彼此显著不同的构象。让我们试试丙烷。

from rdkit import Chem

generator = dc.utils.ConformerGenerator(max_conformers=5)
propane_mol = generator.generate_conformers(Chem.MolFromSmiles('CCC'))
print("Number of available conformers for propane: ", len(propane_mol.GetConformers()))

它只找到了一个构象体。这并不奇怪，因为丙烷是一种非常小的分子，几乎没有任何灵活性。我们尝试再加一个碳。

butane_mol = generator.generate_conformers(Chem.MolFromSmiles('CCCC'))
print("Number of available conformers for butane: ", len(butane_mol.GetConformers()))

现在我们可以为分子创建库仑矩阵。

coulomb_mat = dc.feat.CoulombMatrix(max_atoms=20)
features = coulomb_mat(propane_mol)
print(features)

注意，许多元素都是0。为了将多个分子在一个批量中结合在一起，我们需要把所有的库仑矩阵都设为相同的大小，即使分子的原子数不同。我们指定了 max_atoms=20，因此返回的矩阵的大小为(20,20)。分子只有11个原子，所以只有11 × 11的子矩阵是非零的。

CoulombMatrixEig#

库仑矩阵的一个重要特征是它们不受分子旋转和平动的影响，因为原子间的距离和原子序数不改变。像这样尊重对称性使学习更容易。旋转一个分子并不改变它的物理性质。如果特征化确实发生了变化，那么模型将被迫认识到旋转并不重要，但如果特征化是不变的，那么模型将自动获得该属性。

库仑矩阵在另一个重要的对称性下是变的：原子指标的排列。分子的物理性质与我们称之为“原子1”的原子无关，但库仑矩阵与之相关。为了解决这个问题，引入了 CoulumbMatrixEig 特征器，它使用库仑矩阵的特征值谱，对原子指标的随机排列是不变的。这种特征化的缺点是它包含的信息少得多(N 特征值而不是 N×N 矩阵)，因此模型可以学习的内容将受到更多限制。

CoulombMatrixEig 继承 CoulombMatrix ，通过首先计算分子的不同构象的库仑矩阵，然后计算每个库仑矩阵的特征值来表征一个分子。然后这些特征值被填补以适应各个分子中原子数量的变化。

coulomb_mat_eig = dc.feat.CoulombMatrixEig(max_atoms=20)
features = coulomb_mat_eig(propane_mol)
print(features)

使用拆分器（Splitters）#

在使用机器学习时，通常将数据分为训练集、验证集和测试集。MoleculeNet 加载器会自动完成这个操作。但是应该如何拆分数据呢？这个问题看似简单，但其实很复杂。拆分数据的方法有很多种，选择哪一种方法对结果的可靠性有很大影响。本教程介绍 DeepChem 提供的一些拆分方法。

拆分器（Splitters）#

在 DeepChem 中，将样本拆分为多个数据集的方法是由 Splitter 对象定义的。为数据选择合适的方法非常重要。否则，你训练过的模型可能看起来比实际工作得更好。

想想一个典型的药物开发流程。你可能从筛选成千上万的分子开始，测试它们是否与你感兴趣的目标结合。一旦你找到了一个似乎有效的，你就会尝试通过测试它的成千上万个变体来优化它，寻找一个具有更强亲合力的分子。再然后你可能会在动物身上测试，看它会不会有不可接受的毒性，所以你会尝试更多的变体来解决问题。

这对化学数据集有一个重要的影响：它们通常包含大量彼此非常相似的分子。如果以一种简单的方式将数据分割为训练集和测试集，训练集将包含许多与测试集中的分子非常相似的分子，即使它们并不完全相同。因此，该模型可能在测试集上表现得非常好，但当您试图将其用于与训练数据不太相似的其他数据时，就会失效。

让我们看看在 DeepChem 中的一些分离器。

RandomSplitter#

这是最简单的拆分器之一。它只是以完全随机的方式为训练集、验证集和测试集选择样本。

我们不是说过这是个坏主意吗？这取决于你的数据。如果每个样本都是真正独立于其他样本的，那么这是一个拆分数据的好方法。拆分器没有通用的最佳选择。这完全取决于你的数据集，对于某些数据集，这是一个很好的选择。

RandomStratifiedSplitter#

一些数据集非常不平衡：所有样本中只有一小部分是正样本。在这种情况下，随机拆分有时可能导致验证集或测试集对某些任务只有很少甚至没有正样本。这使得我们无法评估性能。

RandomStratifiedSplitter 通过将正样本和负样本平均拆分来解决这个问题。如果你要求80/10/10的拆分，那么验证集和测试集将不仅包含10%的样本，还包含10%的正样本。

ScaffoldSplitter#

这个拆分器试图解决上面讨论的许多分子彼此非常相似的问题。它识别构成每个分子核心的骨架，并确保具有相同骨架的所有分子都被放入相同的数据集。这仍然不是一个完美的解决方案，因为两个分子可能具有不同的骨架，但在其他方面非常相似，但通常比随机拆分有很大的改进。

ButinaSplitter#

这是另一种尝试解决相似分子问题的拆分器。它根据它们的分子指纹对它们进行聚类，因此具有相似指纹的往往会在同一个数据集中。这种拆分算法所需的时间是跟分子数量的平方成正比，因此它主要用于中小型数据集。

SpecifiedSplitter#

这个拆分器把一切都留给用户。你需要确切地告诉它在每个数据集中放入哪些样本。如果你预先知道特定的拆分适合于你的数据，那么这将非常有用。

一个例子是根据时间拆分。假设有一个研究项目，你再不断地生成和测试新分子。当你获得更多的数据时，你定期地在稳定增长的数据集上重新训练你的模型，然后用它来预测其他尚未测试的分子的结果。验证这种方法是否有效的一种好方法是选择一个特定的截止日期，在当时拥有的所有数据上训练模型，并看看它对以后生成的其他数据的预测效果如何。

使用不同拆分器的效果#

让我们来看一个例子。我们将使用 random、scaffold、和 Butina 拆分器加载Tox21毒性数据集。对于每一个，我们训练一个模型，并在训练集和测试集上对其进行评估。

import deepchem as dc

splitters = ['random', 'scaffold', 'butina']
metric = dc.metrics.Metric(dc.metrics.roc_auc_score)
for splitter in splitters:
    tasks, datasets, transformers = dc.molnet.load_tox21(featurizer='ECFP', split=splitter)
    train_dataset, valid_dataset, test_dataset = datasets
    model = dc.models.MultitaskClassifier(n_tasks=len(tasks), n_features=1024, layer_sizes=[1000])
    model.fit(train_dataset, nb_epoch=10)
    print('splitter:', splitter)
    print('training set score:', model.evaluate(train_dataset, [metric], transformers))
    print('test set score:', model.evaluate(test_dataset, [metric], transformers))
    print()

它们在训练集上产生非常相似的性能，但 random 拆分器在测试集上有更高的性能。Scaffold 拆分器有较低的测试集得分，Butina 拆分器甚至更低。这是否意味着 random 拆分器更好？不!这意味着 random 拆分器并不能准确地衡量模型的工作情况。因为测试集包含许多与训练集中的分子非常相似的分子，所以它不是真正独立的。它使模型看起来比实际工作得更好。Scaffold 和 Butina 拆分器可以更好地说明未来独立数据的预期结果。

高级模型训练#

在到目前为止的教程中，我们遵循了一个简单的训练模型的过程：加载数据集，创建模型，调用 fit() ，对其进行评估，并认为我们完成了。对于一个例子来说，这很好，但在真正的机器学习项目中，这个过程通常更复杂。在本教程中，我们将讨论训练模型的更实际的工作流程。

超参数优化#

让我们从加载 HIV 数据集开始。它根据是否抑制艾滋病毒复制对超过4万个分子进行了分类。

import deepchem as dc

tasks, datasets, transformers = dc.molnet.load_hiv(featurizer='ECFP', split='scaffold')
train_dataset, valid_dataset, test_dataset = datasets

现在让我们用它来训练一个模型。我们将使用 MultitaskClassifier ，它只是一个密集层的堆栈。但是仍然有很多的选择。应该有多少层，每层应该有多少个神经元？我们应该使用什么 dropout rate？学习率是多少？

这些被称为超参数。选择它们的标准方法是大量地尝试，在训练集中训练每个模型，并在验证集中评估它。这让我们知道哪种方法效果最好。

你可以用手来做，但通常让电脑帮你做更容易。DeepChem 提供了一系列超参数优化算法，这些算法在 dc.hyper 库里。对于这个例子，我们将使用 GridHyperparamOpt ，这是最基本的方法。我们只是给它一个超参数选项的列表，它就会一个一个地尝试它们的所有组合。

选项列表由我们提供的 dict 定义。对于模型的每个参数，我们提供一个列表来尝试。在这个例子中，我们考虑了三组可能的隐藏层：神经元为500的单层，神经元为1000的单层，或者神经元为1000的两层。我们还考虑了两个 drop out rate (20%和50%)和两个学习率(0.001和0.0001)。

params_dict = {
    'n_tasks': [len(tasks)],
    'n_features': [1024],
    'layer_sizes': [[500], [1000], [1000, 1000]],
    'dropouts': [0.2, 0.5],
    'learning_rate': [0.001, 0.0001]
}
optimizer = dc.hyper.GridHyperparamOpt(dc.models.MultitaskClassifier)
metric = dc.metrics.Metric(dc.metrics.roc_auc_score)
best_model, best_hyperparams, all_results = optimizer.hyperparam_search(
        params_dict, train_dataset, valid_dataset, metric, transformers)

hyperparam_search() 返回三个参数:它找到的最佳模型、该模型的超参数，以及每个模型的验证得分的完整列表。让我们来看看最后一个。

print(all_results)

我们可以看到一些通用规律。使用两个学习率较高的层并不十分有效。似乎更深的模型需要更小的学习率。我们还发现，20%的 dropout 通常比50%好。一旦我们根据这些观察结果缩小模型列表，所有验证分数都非常接近，可能接近到足以使剩余的变化主要是噪声。我们使用的剩余超参数集似乎没有太大区别，所以让我们任意选择一个具有神经元数为1000的单层，学习率为0.0001的模型。

提前停止#

还有一个我们尚未考虑的重要超参数：我们训练模型的时间 GridHyperparamOpt 对每个模型进行固定的、相当少的轮次（epochs）训练。这不一定是最好的数字。

你可能认为你训练的时间越长，你的模型就会越好，但这通常不是真的。如果训练时间过长，模型通常会过度拟合训练集的无关细节。你可以判断何时发生这种情况，因为验证集分数停止增加，甚至可能会减少，而训练集的分数继续提高。

幸运的是，我们不需要为不同的步数训练许多次模型来确定最佳步数。我们只需要训练它一次，监视验证分数，并保留使它最大化的任何参数。这就是所谓的“提前停止”。DeepChem 的 ValidationCallback 类可以自动为我们完成这个操作。在下面的例子中，我们要求它每1000个训练步骤计算验证集的ROC AUC。如果你添加 save_dir 参数，它也会将最佳模型参数保存到硬盘。

model = dc.models.MultitaskClassifier(n_tasks=len(tasks),
                                    n_features=1024,
                                    layer_sizes=[1000],
                                    dropouts=0.2,
                                    learning_rate=0.0001)
callback = dc.models.ValidationCallback(valid_dataset, 1000, metric)
model.fit(train_dataset, nb_epoch=50, callbacks=callback)

学习率计划(Learning Rate Schedules)#

在上面的例子中，我们在整个训练过程中使用了固定的学习率。在某些情况下，在训练过程中改变学习速度的效果更好。要在 DeepChem 中做到这一点，我们只需为 learning_rate 参数指定一个 LearningRateSchedule 对象而不是一个数字。在下面的例子中，我们会使用一个会指数下降的学习率。它从0.0002开始，然后在每走1000步之后乘以0.9。

learning_rate = dc.models.optimizers.ExponentialDecay(0.0002, 0.9, 1000)
model = dc.models.MultitaskClassifier(n_tasks=len(tasks),
                                    n_features=1024,
                                    layer_sizes=[1000],
                                    dropouts=0.2,
                                    learning_rate=learning_rate)
model.fit(train_dataset, nb_epoch=50, callbacks=callback)

从实验数据创建一个高精确度的模型#

在本教程中，我们将了解从实验数据创建一个新的Dataset所涉及的内容。我们将看到，创建Dataset对象的机制只是整个过程的一小部分。大多数真实的数据集在适合用于训练模型之前都需要大量的清理和质量检查。

使用数据文件#

假设你得到了一个你的实验合作伙伴收集的数据。你可能想用这些数据构建一个机器学习模型。

你如何将这些数据转换为能够创建有用模型的数据集？

从新数据构建模型可能会带来一些挑战。也许数据记录的方式不方便。此外，数据可能包含噪声。例如，由于大量的外部变量以及与收集多个样本相关的难度和成本，这在生物检测中是常见的。这是一个问题，因为你不希望你的模型拟合于这种噪声。

因此，主要有两个挑战：

解析数据
数据去噪

在本教程中，我们将通过一个例子，使用一个有药物实验数据的 excel 电子表格创建一个数据集。在我们深入这个例子之前，让我们简单回顾一下 DeepChem 的输入文件处理和特征化功能。

输入格式#

DeepChem 支持各种输入文件。例如，可接受的文件格式包括 .csv、.sdf、.fasta、.png、.tif 和其他文件格式。特定文件格式的加载由与该格式关联的 Loader 类控制。例如，要加载一个 .csv 文件，我们使用 CSVLoader 类。下面是一个符合 CSVLoader 要求的 .csv文件示例。

包含SMILES字符串的列。
包含实验测量值的列。
(可选)包含唯一化合物标识符的列。

下面是一个输入文件示例。

化合物ID	log（溶解度）测量值，单位为摩尔/升	smiles
苯并噻唑	-1.5	c2ccc1scnc1c2

这里的“smiles”列包含 SMILES 字符串，“log（溶解度）测量值，单位为摩尔/升”包含实验测量值，“化合物ID”包含唯一的化合物标识符。

数据特征化#

大多数机器学习算法要求输入数据为向量。然而，药物发现数据集的输入数据通常以分子列表和相关实验数据的形式出现。为了加载数据，我们使用 dc.data.DataLoader 的一个子类。例如 dc.data.CSVLoader 或 dc.data.SDFLoader 。用户可以子类化 dc.data.DataLoader 加载任意文件格式。所有的加载器都要求 dc.feat.Featurizer 对象作为参数，它指定如何将分子转换为向量。DeepChem 提供了 dc.feat.Featurizer 的许多不同子类。

解析数据#

为了读取数据，我们将使用pandas数据分析库。

为将药物名称转换为smiles字符串，我们将使用 pubchempy。这不是一个标准的 DeepChem 依赖，但你可以用 conda install pubchemy 来安装这个库。

import os
import pandas as pd
from pubchempy import get_cids, get_compounds

Pandas 很神奇，但它不能自动知道在哪里找到你感兴趣的数据。你可能必须首先使用GUI查看数据集。

现在，我们来看看由 LibreOffice 渲染的这个数据集的截图。

为此，我们将导入 Image 和 os。

import os
from IPython.display import Image, display
current_dir = os.path.dirname(os.path.realpath('__file__'))
data_screenshot = os.path.join(current_dir, 'assets/dataset_preparation_gui.png')
display(Image(filename=data_screenshot))

https://github.com/deepchem/deepchem/raw/master/examples/tutorials/assets/dataset_preparation_gui.png

我们可以在第二页看到感兴趣的数据，包含在“TA ID”、“N #1(%)”和“N #2 (%)”列中。

此外，该电子表格的大部分格式似乎是为了便于人类阅读(多列标题，带有空格和符号的列标签，等等)。这使得创建一个整洁的数据框架对象变得更加困难。因此，我们将砍掉一切不必要的或不方便的东西。

raw_data_file = '..\..\datasets\Positive Modulators Summary_ 918.TUC _ v1.xlsx'
raw_data_excel = pd.ExcelFile(raw_data_file)

# second sheet only
raw_data = raw_data_excel.parse(raw_data_excel.sheet_names[1])

# preview 5 rows of raw dataframe
print(raw_data.loc[raw_data.index[:5]])

注意，实际的行标题存储在第1行，而不是上面的第0行。

# remove column labels (rows 0 and 1), as we will replace them
# only take data given in columns "TA ID" "N #1 (%)" (3) and "N #2 (%)" (4)
raw_data = raw_data.iloc[2:, [2, 6, 7]]

# reset the index so we keep the label but number from 0 again
raw_data.reset_index(inplace=True)

## rename columns
raw_data.columns = ['label', 'drug', 'n1', 'n2']

# preview cleaner dataframe
raw_data

这种格式更接近我们所需要的。

现在，让我们用药物名称并获取它们的 smiles 字符串(DeepChem 所需的格式)。

drugs = raw_data['drug'].values

对于其中的大多数，我们可以通过 get_compounds 对象(使用 pubchempy)的 canonical_smiles 属性获得smiles字符串。

get_compounds(drugs[1], 'name')
get_compounds(drugs[1], 'name')[0].canonical_smiles

然而，其中一些药物名称有不同大小的空间和符号(·、(±)等)，pubchempy 可能无法读取这些名称。

对于这个任务，我们将通过普通表达式进行一些修改。同时，我们注意到所有离子都是以缩写形式写的，因此需要展开。因此，我们使用字典，将缩短的离子名称映射到 pubchempy 可识别的版本。

不幸的是，你可能会遇到一些极端情况，需要进行更多的修改。

import re

ion_replacements = {
    'HBr': ' hydrobromide',
    '2Br': ' dibromide',
    'Br': ' bromide',
    'HCl': ' hydrochloride',
    '2H2O': ' dihydrate',
    'H20': ' hydrate',
    'Na': ' sodium'
}

ion_keys = ['H20', 'HBr', 'HCl', '2Br', '2H2O', 'Br', 'Na']

def compound_to_smiles(cmpd):
    # remove spaces and irregular characters
    compound = re.sub(r'([^\s\w]|_)+', '', cmpd)

    # replace ion names if needed
    for ion in ion_keys:
        if ion in compound:
            compound = compound.replace(ion, ion_replacements[ion])

    # query for cid first in order to avoid timeouterror
    cid = get_cids(compound, 'name')[0]
    smiles = get_compounds(cid)[0].canonical_smiles

    return smiles

现在让我们把所有这些化合物转化为 smiles。这个转换将需要几分钟，所以这可能是一个去喝杯咖啡或茶，并在运行的时候休息的不错的时候！请注意，这种转换有时会失败，因此我们在下面添加了一些错误处理来捕获这种情况。

smiles_map = {}
for i, compound in enumerate(drugs):
    try:
        smiles_map[compound] = compound_to_smiles(compound)
    except:
        print("Errored on %s" % i)
        continue

smiles_data = raw_data
# map drug name to smiles string
smiles_data['drug'] = smiles_data['drug'].apply(lambda x: smiles_map[x] if x in smiles_map else None)

# preview smiles data
smiles_data.loc[smiles_data.index[:5]]

Hooray，我们已经将每个药物名称映射到对应的 smiles 代码。

现在，我们需要查看数据并尽可能多地去除噪声。

数据去噪（De-noising data）#

在机器学习中，我们知道没有免费的午餐。你需要花时间分析和理解你的数据，以便构建你的问题并确定合适的模型框架。从这个过程中收集到的结论决定数据的处理方式。

要问自己的问题：

你想要完成什么？
你的化验结果是什么？
数据的结构是什么？
数据是否有意义？
以前尝试过什么？

对于本项目：

我想建立一个能够预测任意小分子药物对特定离子通道蛋白亲和力的模型
对输入的某个药物，数据描述通道抑制能力
几百种 n = 2 的药物
需要更仔细地观察数据集*
对于这种蛋白质，没有信息
这将涉及绘图，所以我们将导入 matplotlib 和 seaborn。我们还需要研究分子结构，因此我们将导入 rdkit。我们还将使用 seaborn 库，你可以使用 conda install seaborn 来安装它。

import matplotlib.pyplot as plt
%matplotlib inline

import seaborn as sns
sns.set_style('white')

from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem import Draw, PyMol, rdFMCS
from rdkit.Chem.Draw import IPythonConsole
from rdkit import rdBase
import numpy as np

我们的目标是建立一个小分子模型，所以我们要确保我们的分子都很小。这可以用每个 smiles 字符串的长度来近似表示。

smiles_data['len'] = [len(i) if i is not None else 0 for i in smiles_data['drug']]
smiles_lens = [len(i) if i is not None else 0 for i in smiles_data['drug']]
sns.histplot(smiles_lens)
plt.xlabel('len(smiles)')
plt.ylabel('probability')

有些看起来相当大， len(smiles) > 150。让我们看看它们长什么样。

# indices of large looking molecules
suspiciously_large = np.where(np.array(smiles_lens) > 150)[0]

# corresponding smiles string
long_smiles = smiles_data.loc[smiles_data.index[suspiciously_large]]['drug'].values

# look
Draw._MolsToGridImage([Chem.MolFromSmiles(i) for i in long_smiles], molsPerRow=1)

正如怀疑的那样，这些不是小分子，所以我们将从数据集中删除它们。这里的论点是，这些分子可以作为抑制剂仅仅是因为它们很大。它们更有可能在空间上阻塞通道，而不是扩散到内部并结合(这是我们感兴趣的)。

这里的经验是删除不适合你的用例的数据。

# drop large molecules
smiles_data = smiles_data[~smiles_data['drug'].isin(long_smiles)]

现在，让我们看看数据集的数值结构。

首先，检查 NaN。

nan_rows = smiles_data[smiles_data.isnull().T.any().T]
print(nan_rows[['n1', 'n2']])

我不相信n=1，所以我不要它们了。

然后，让我们检验n1和n2的分布。

df = smiles_data.dropna(axis=0, how='any')
# seaborn jointplot will allow us to compare n1 and n2, and plot each marginal
sns.jointplot(x='n1', y='n2', data=smiles_data)

我们看到大部分数据都包含在高斯团中，中心略低于零。我们可以看到，在左下方有几个明显有活性的数据点，右上方有一个。这些都与大多数数据有所区别。我们如何处理团中的数据?

因为n1和n2表示相同的测量值，理想情况下它们的值应该是相同的。这个图应该与对角线紧密对齐，皮尔逊相关系数应该是1。我们看到事实并非如此。这有助于我们对实验的误差有一个概念。

让我们更仔细地观察误差，绘制(n1-n2)的分布。

diff_df = df['n1'] - df['n2']

sns.histplot(diff_df)
plt.xlabel('difference in n')
plt.ylabel('probability')

这看起来很像高斯分布，让我们通过 scipy 拟合高斯分布得到95%的置信区间，并取2个标准差。

from scipy import stats
mean, std = stats.norm.fit(np.asarray(diff_df, dtype=np.float32))
ci_95 = std*2
print(ci_95)

现在，我不相信置信区间之外的数据，因此会从df中删除这些数据点。

例如，在上图中，至少有一个数据点的n1-n2 > 60。这是令人不安的。

noisy = diff_df[abs(diff_df) > ci_95]
df = df.drop(noisy.index)
sns.jointplot(x='n1', y='n2', data=df)

现在数据看起来好多了!

我们求n1和n2的平均值，取误差条(error bar)为ci_95。

avg_df = df[['label', 'drug']].copy()
n_avg = df[['n1', 'n2']].mean(axis=1)
avg_df['n'] = n_avg
avg_df.sort_values('n', inplace=True)

现在，让我们看一下带有误差条的排序数据。

plt.errorbar(np.arange(avg_df.shape[0]), avg_df['n'], yerr=ci_95, fmt='o')
plt.xlabel('drug, sorted')
plt.ylabel('activity')

现在，我们来鉴定活性化合物。

在我的情况下，这需要领域知识。因为我在这一领域工作过，也咨询过这方面的专家，所以我对活性绝对值大于25的化合物很感兴趣。这与我们想要建模的期望药效有关。

如果你不确定如何在活动和非活动之间划清界限，则此界限可能被视为超参数。

actives = avg_df[abs(avg_df['n'])-ci_95 > 25]['n']

plt.errorbar(np.arange(actives.shape[0]), actives, yerr=ci_95, fmt='o'

# summary
print(raw_data.shape, avg_df.shape, len(actives.index))

综上所述，我们有：

删除了与我们希望回答的问题无关的数据（仅保留小分子数据）
删除了NaNs
确定了我们测量值的噪声
删了除异常嘈杂的数据点
确定了活性物（使用领域知识来确定一个阈值）

确定模型类型、数据集的最终形式和健全加载（sanity load）#

现在，我们应该使用什么模型框架呢?

鉴于我们有392个数据点和6个活性物,这些数据将被用来建立一个被称为（ low data one-shot classifier）的分类器(10.1021/acscentsci.6b00367)。如果有类似性质的数据集，就可能使用迁移学习，但目前情况并非如此。

让我们将逻辑应用到数据中，以便将其转换为适合分类的二进制格式。

# 1 if condition for active is met, 0 otherwise
avg_df.loc[:, 'active'] = (abs(avg_df['n'])-ci_95 > 25).astype(int)

现在，将其保存到文件中。

avg_df.to_csv('modulators.csv', index=False)

现在，我们将把这个数据框架转换为一个 DeepChem 数据集。

import deepchem as dc
dataset_file = 'modulators.csv'
task = ['active']
featurizer_func = dc.feat.ConvMolFeaturizer()

loader = dc.data.CSVLoader(tasks=task, feature_field='drug', featurizer=featurizer_func)
dataset = loader.create_dataset(dataset_file)

最后，以某种方式对数据进行数值变换通常是有利的。例如，有时对数据进行归一化或将均值归零是很有用的。这取决于手头的任务。

DeepChem 内置了许多有用的变换器（transformer），位于 deepchem.transformers.transformers 类库。

因为这是一个分类模型，而且活性药物的数量很低，所以我将使用一个平衡变换器。当我开始训练模型时，我把这个变换器当作一个超参数。它被证明无疑地提高了模型的性能。

transformer = dc.trans.BalancingTransformer(dataset=dataset)
dataset = transformer.transform(dataset)

现在，让我们将平衡的数据集对象保存到硬盘，然后重新加载它作为健全性检查。

dc.utils.save_to_disk(dataset, 'balanced_dataset.joblib')
balanced_dataset = dc.utils.load_from_disk('balanced_dataset.joblib')

完。

把多任务机器学习模型应用在工作中#

[1] https://pubs.acs.org/doi/10.1021/ci8002649

把多任务机器学习模型应用在工作中#

本笔记本介绍了在MUV[1]上创建多任务模型的过程。目的是演示多任务方法如何在数据很少或非常不平衡的情况下提高性能。

MUV 数据集是一个具有挑战性的分子设计基准，它由17个不同的“靶标”组成，每个靶标只有几个“活性”化合物。一共有93087种化合物，但没有一个任务的活性化合物超过30种，许多任务的活性化合物甚至更少。用如此少的正样本训练一个模型是非常具有挑战性的。多任务模型通过训练一个单一模型来解决这个问题，该模型可以同时预测所有不同的目标。如果一个特性对预测一个任务有用，那么它通常也对预测其他几个任务有用。每增加一个任务，就可以更容易地学习重要的特性，从而提高其他任务的性能[2]。

首先，让我们加载 MUV 数据集。MoleculeNet 加载器函数自动将其分为训练集、验证集和测试集。由于正样本太少，我们使用 stratified 拆分器来确保测试集有足够的正样本进行评估。

import deepchem as dc
import numpy as np

tasks, datasets, transformers = dc.molnet.load_muv(split='stratified')
train_dataset, valid_dataset, test_dataset = datasets

现在让我们用它来训练一个模型。我们将使用 MultitaskClassifier，它是一个完全连接层的简单堆栈。

n_tasks = len(tasks)
n_features = train_dataset.get_data_shape()[0]
model = dc.models.MultitaskClassifier(n_tasks, n_features)
model.fit(train_dataset)

让我们看看它在测试集上的表现如何。我们循环这17个任务，并为每个任务计算 ROC AUC。

y_true = test_dataset.y
y_pred = model.predict(test_dataset)
metric = dc.metrics.roc_auc_score
for i in range(n_tasks):
    score = metric(dc.metrics.to_one_hot(y_true[:,i]), y_pred[:,i])
    print(tasks[i], score)

还不错!回想一下，随机猜测会产生0.5的ROC AUC分数，而一个完美的预测器会得到1.0分。大多数任务比随机猜测的表现要好得多，其中许多任务的得分都在0.9以上。

参考书目#

[2] https://pubs.acs.org/doi/abs/10.1021/acs.jcim.7b00146

建立蛋白质-配体相互作用模型#

在本教程中，我们将带你了解如何使用机器学习和分子对接方法来预测蛋白质配体复合体的结合能。回想一下，配体是一种与蛋白质相互作用的小分子(通常是非共价的)。分子对接通过几何计算找到一种“结合构象”，小分子与蛋白质在合适的结合口袋中相互作用(也就是蛋白质上有凹槽的区域，小分子可以进入其中)。

蛋白质的结构可以通过像 Cryo-EM 或 X 射线晶体学这样的技术来实验鉴定。这可以是用于基于结构的药物发现的一个强大工具。关于对接的更多内容，请阅读 AutoDock Vina paper 和 deepchem.dock 文档。有许多图形用户界面和命令行界面程序(如 AutoDock)用于执行分子对接。在这里，我们要展示如何通过 DeepChem 编程实现对接，从而能使自动化，和轻松与机器学习流程整合。

通过本教程,你将沿着包括下列的事项学习：

加载蛋白质-配体复合物数据集 [PDBbind]
通过编程执行分子对接
使用相互作用指纹特征化蛋白质-配体复合物
拟合随机森林模型并预测结合亲和力

为了开始本教程，我们将使用一个简单的预处理数据集文件，它以 gzip 文件的形式出现。每一行都是一个分子系统，每一列都代表了关于这个系统的不同信息。例如，在这个例子中，每一行都反映一个蛋白质-配体复合物，然后每一个列是：复合物唯一标识符；配体的 SMILES 字符串；配体与复合物中蛋白质的结合亲和力 (Ki)；所有行的只有蛋白质的PDB文件 Python 列表；以及所有行的只有配体的PDB文件 Python 列表。

准备#

蛋白质-配体复合物数据#

在进行对接时，蛋白质和配体的可视化是很有帮助的。不幸的是，谷歌Colab（如果要使用谷歌Colab，国内需要翻墙）目前不支持我们进行可视化所需要的 Jupyter 小部件。在本地机器上可以安装 MDTraj 和 nglview 来查看我们正在使用的蛋白质-配体复合物。

pip install -q mdtraj nglview
# jupyter-nbextension enable nglview --py --sys-prefix  # for jupyter notebook
# jupyter labextension install  nglview-js-widgets  # for jupyter lab


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

为了演示对接过程，这里我们将使用一个 csv，其中包含来自 PDBbind 的配体的 SMILES 字符串以及配体和蛋白质靶点的 PDB 文件。稍后，我们会使用标签列训练一个模型来预测结合亲和力。我们还将从零开始展示如何下载和特征化 PDBbind 以训练模型。

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']]

让我们看看 raw_dataset 是什么样的：

raw_dataset.head(2)

修复 PDB 文件#

接下来，让我们获取一些用于可视化和对接的 PDB 蛋白质文件。我们将使用来自 raw_dataset 的 PDB ID，并使用 pdbfixer 直接从蛋白质数据库下载 PDB 文件。我们还将使用 RDKit 对结构进行清理。这确保蛋白质和配体文件中存在的任何问题（非标准残基，化学有效性等）会得到纠正。请随意修改这些框和 pbid，以讨论新的蛋白质-配体复合物。我们在这里注意到 PDB 文件是复杂的，需要人类的判断来准备对接的蛋白质结构。DeepChem 包含许多对接程序来帮助你准备蛋白质文件，但在尝试对接之前应该检查准备结果。

from openmm.app import PDBFile
from pdbfixer import PDBFixer

from deepchem.utils.vina_utils import prepare_inputs

# consider one protein-ligand complex for visualization
pdbid = raw_dataset['pdb_id'].iloc[1]
ligand = raw_dataset['smiles'].iloc[1]

%%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))

可视化#

如果你在 Colab 之外，你可以执行下面的代码，并使用 MDTraj 和 MDTraj 来可视化蛋白质和配体。

import mdtraj as md
import nglview

from IPython.display import display, Image

让我们来看看数据集中的第一个蛋白质-配体对：

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

我们将使用函数 nglview. show_mdtraj 来查看我们的蛋白质和配体。注意，只有当你安装了nglview并启用必要的笔记本扩展时，这才会起作用。

v = nglview.show_mdtraj(ligand_mdtraj)

display(v)  # interactive view outside Colab

现在我们已经知道了配体的样子，让我们看看我们的蛋白质：

view = nglview.show_mdtraj(protein_mdtraj)
display(view)  # interactive view outside Colab

分子对接#

好了，现在我们已经有了数据和基本的可视化工具，让我们看看是否可以使用分子对接来估计蛋白质配体系统之间的结合亲和力。

设置对接任务有三个步骤，你应该尝试不同的设置。我们需要明确的三件事是：

如何识别目标蛋白质中的结合口袋；
如何生成结合口袋中配体的取向(几何构象);
如何“评分”一个构象。

记住，我们的目标是识别与目标蛋白强烈相互作用的候选配体，这可以通过评价分数反映出来。

DeepChem 有一种简单的内置方法，可以识别蛋白质中的结合口袋。它是基于凸面外壳法（convex hull method ）的。该方法的工作原理是在蛋白质结构周围创建一个三维多面体（convex hull），并确定最接近凸面外壳的蛋白质表面原子。由于考虑了一些生物化学性质，所以该方法不是纯几何的。它的优点是计算成本低，足以满足我们的目的。

finder = dc.dock.binding_pocket.ConvexHullPocketFinder()
pockets = finder.find_pockets('3cyx.pdb')
len(pockets)  # number of identified pockets

构象生成相当复杂。幸运的是，使用 DeepChem 的基于 AutoDock Vina 引擎的构象生成器使我们能够快速启动和运行构象生成。

vpg = dc.dock.pose_generation.VinaPoseGenerator()

我们可以从 deepchem.dock. pose_scoring 中指定一个包括排斥和疏水相互作用和氢键的构象评分函数。Vina 将帮我们处理处理这个问题，所以我们将允许 Vina 为构象计算分数。

mkdir -p vina_test
%%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 将以 kcal/mol 为单位返回9个能量最低的构象。

print(scores)

我们能同时观察蛋白质和配体的复合物吗？是的，但我们需要把这些分子组合成一个 RDkit 分子。

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

现在我们来显现一下复合体。我们可以看到配体插入到蛋白质的一个口袋里。

v = nglview.show_rdkit(complex_mol)
display(v)

现在我们已经了解了整个过程的各个部分，我们可以使用 DeepChem 的 Docker 类将它们组合在一起。Docker 将创建一个生成器，生成复合结构和对接分数组成的元组。

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 中的“精制（refined）”集获取4852个蛋白质-配体复合物或获取整个“一般（general）”集。为了简单起见，我们将坚持使用我们已经处理过的大约100个复合物来训练我们的模型。

接下来，我们需要一种方法，将我们的蛋白质-配体复合物转换成可以被学习算法使用的表示形式。理想情况下，我们应该有神经网络蛋白-配体复合体指纹，但 DeepChem 还没有这种良好的机器学习指纹。然而，我们确实有手动调整好的特征器，可以帮助我们在这里的挑战。

在接下来的教程中，我们将使用两种类型的指纹， CircularFingerprint 和 ContactCircularFingerprint 。DeepChem 还拥有体素化器（voxelizers）和网格描述符（grid descriptors），可将包含原子排列的 3D 体块转换为指纹。这些特征器对于理解蛋白质-配体复合物非常有用，因为它们允许我们将复合物转换为可以传递到简单机器学习算法中的向量。首先，我们要创建 CircularFingerprints 。它们将小分子转化为片段向量。

pdbids = raw_dataset['pdb_id'].values
ligand_smiles = raw_dataset['smiles'].values


%%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))

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 ID，并删除任何没有相应配体的蛋白质。

# 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')

len(proteins), len(ligands)

pdbids = [f[:-4] for f in proteins]
small_dataset = raw_dataset[raw_dataset['pdb_id'].isin(pdbids)]
labels = small_dataset.label

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

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

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 的精致集，请取消注释并享受一杯咖啡。否则，你可以继续使用我们上面构造的小数据集。

# # 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) task_type， (b) model_params，另一个 `dict` 如下所示，以及 (c) 一个 `model_instance` 定义你想要的模型类型，在本例中是 `RandomForestRegressor` 。

from sklearn.ensemble import RandomForestRegressor

from deepchem.utils.evaluate import Evaluator
import pandas as pd

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)

注意，测试集的 \(R^2\) 值表明模型没有产生有意义的输出。事实证明，预测结合亲和力是 困难的 。本教程并不是要展示如何创建最先进的预测结合亲和力的模型，而是为你提供使用分子对接、特征化复合物和训练模型生成自己的数据集的工具。

# 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"]))

我们使用的是非常小的数据集和过于简单的表示，所以测试集的性能非常糟糕也就不足为奇了。

# Compare predicted and true values
list(zip(model.predict(train_dataset), train_dataset.y))[:5]

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

蛋白质-配体复合物的显现#

在上一节中，我们只特征化了配体。这一次，让我们看看能否利用我们的结构信息，对蛋白质-配体指纹做些有意义的事情。首先，我们需要重新特征化数据集，但这次使用接触圆形指纹（contact circular fingerprint）。

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)

现在让我们在这个数据集上训练一个简单的随机森林模型。

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)

让我们看看我们的准确性是什么样的！

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"]))

好的，看起来我们的精度比仅有配体数据集要低。尽管如此，拥有一个蛋白质-配体模型可能还是有用的，因为它可能比纯配体模型学习到不同的特征。

使用原子卷积网络建立蛋白质配体相互作用模型#

本 DeepChem 教程要介绍原子卷积神经网络（Atomic Convolutional Neural Network, ACNN）。我们将看到 AtomicConvModel 的架构，并编写一个简单的程序来运行原子卷积（Atomic Convolutions）。

背景知识：#

在数学中，一个距离矩阵是一个包含一组点两两之间距离的矩阵(即二维数组)。因此给定N个欧几里得空间中的点，其距离矩阵就是一个非负实数作为元素的N×N的对称矩阵。这些点两两之间点对的数量，N×(N-1)/2，也就是距离矩阵中独立元素的数量。—— 百度百科

邻近列表程序（neighbor list routine）：基本思路是，把一个分子放在一个长方体容器，然后容器被分割成许多小的盒子，最后某个小盒子内的原子的邻近列表由这个小盒子附近的小盒子内的全部原子组成。详细的介绍如下：

首先确定容器的长，宽，高。
然后以固定的长度分割容器，比如：有个正方体容器，边为40埃，每边分成八份，因此总共有8*8*8=512个小盒子，每一个盒子的边长为5埃，最后给每个小盒子制定坐标。比如：把最左边，最下面，最前面的小盒子的坐标设为(1,1,1)，那么、最右边，最上面，最后面的小盒子坐标为(8,8,8)
然后根据每个原子的正交坐标系坐标确定每个小盒子中的原子，我们可以把这些数据保存在一个(x,y,z,i)的四维数列中，其中x,y,z是原子所在的小盒子的坐标（如上面所述），i为原子编号。
上面所述的四维数列，最开始以源自编号排序的，后面根据原子所在的小盒子坐标(x,y,z)排升序。
人为的给定某个范围，分子中某个原子的邻近列表是上面所述的排序好的四维数列中围绕着这个原子的所给定范围内的所有原子。

通过上述过程，我们就能得到分子中某个原子的邻近列表。也就是说，某个原子的邻近列表是这个原子所在的小盒子附近的小盒子中的所有原子。当然，“附近”的定义取决于你。我们使用这种程序是因为，这个程序的运行时间短，跟总共原子数乘以M成正比。（这里的M指的是上面所述的范围）。用公式表示的话：O(NM)（N是总原子数）。

参考文献：https://doi.org/10.1002/jcc.540100709

ACNN 架构#

ACNN 直接利用分子的局部三维结构，通过端到端同时优化模型和特征化，分层学习更复杂的化学特征。

原子卷积利用邻近列表距离矩阵从不一定包含空间局部性的输入表示(笛卡尔原子坐标)中提取具有局部化学环境的特征。以下方法用于构建ACNN架构:

距离矩阵#

距离矩阵 R 是由笛卡尔原子坐标 X 构造的。它计算距离张量 D 的距离。距离矩阵构造接受 (N,3) 坐标矩阵 C 作为输入。这个矩阵被“邻近列表”成 (N,M) 矩阵 R。

R = tf.reduce_sum(tf.multiply(D, D), 3)     # D: Distance Tensor
R = tf.sqrt(R)                              # R: Distance Matrix
return R

原子卷积#

原子卷积的输出由距离矩阵 R 和原子序数矩阵 Z 构造。矩阵 R 被送入到步数（stride）为1、深度为 \(N_{at}\) 的 (1x1) 过滤器（filter），其中 \(N_{at}\) 是分子系统中存在的唯一原子序数(原子类型)的数量。原子卷积核是一个作用于邻近距离矩阵R的阶跃函数。

辐射池化层（Radial Pooling layer）#

辐射池化基本上是一个降维过程，对原子卷积的输出进行降采样。降维过程通过特征分类提供一种抽象的表示形式，以及减少需要学习的参数数量，从而防止过拟合。

从数学上讲，辐射池化层在大小为 (1xMx1)、步数为1、深度为 \(N_r\) 的张量切片(接受域)上进行池化，其中 \(N_r\) 是所需辐射过滤器的数量，M是邻近范围。

原子全连接网络（Atomistic fully connected network）#

原子卷积层通过将扁平的辐射池化层输出（\(N, N_{at} \cdot N_r\)）输入到原子卷积进行堆叠。最后，我们将张量按行(按原子)输入到一个完全连接的网络中。对于给定分子中的每个原子，使用相同的全连接权重和偏置。

现在我们已经了解了 ACNN 的架构概述，我们将尝试更深入地研究模型，看看如何训练它，以及我们期望的输出是什么。

对于训练，我们将使用公开可用的 PDBbind 数据集。在这个例子中，每一行表示一个蛋白质-配体复合物，目标结果是配体与复合物中蛋白质的结合亲和力（\(K_i\)）。

import deepchem as dc
import os

import numpy as np
import tensorflow as tf

import matplotlib.pyplot as plt

from rdkit import Chem

from deepchem.molnet import load_pdbbind
from deepchem.models import AtomicConvModel
from deepchem.feat import AtomicConvFeaturizer

获取蛋白质配体数据#

如果你已经完成了关于建立蛋白质-配体相互作用模型的教程 13，你就已经熟悉了如何从 PDBbind 获得一组数据来训练我们的模型。由于我们在上一篇教程中详细探讨了分子复合物，这次我们将简单地初始化一个 AtomicConvFeaturizer 并直接使用 MolNet 加载 PDBbind 数据集。

f1_num_atoms = 100  # maximum number of atoms to consider in the ligand
f2_num_atoms = 1000  # maximum number of atoms to consider in the protein
max_num_neighbors = 12  # maximum number of spatial neighbors for an atom

acf = AtomicConvFeaturizer(frag1_num_atoms=f1_num_atoms,
                    frag2_num_atoms=f2_num_atoms,
                    complex_num_atoms=f1_num_atoms+f2_num_atoms,
                    max_num_neighbors=max_num_neighbors,
                    neighbor_cutoff=4)

load_pdbbind 允许我们指定是要使用整个蛋白质还是只使用结合口袋（ pocket=True ）进行特征化。只使用口袋节省内存和加快特征化。我们还可以使用由约200个高质量复合物组成的“核心（core）”数据集来快速测试我们的模型，或者使用更大的由近5000个复合物组成的的“精炼（refined）”集来实现更多的数据点和更严谨的训练/验证。在Colab上，只需要一分钟就可以特征化核心 PDBbind 数据集！这是非常不可思议的，这意味着你可以快速地试验不同的特征化和模型架构。

%%time
tasks, datasets, transformers = load_pdbbind(featurizer=acf,
                                            save_dir='.',
                                            data_dir='.',
                                            pocket=True,
                                            reload=False,
                                            set_name='core')

不幸的是，如果你试图使用“精炼（refined）”数据集，就会发现有些复合物无法被特征化。要解决这个问题，除了增加 complex_num_atoms ，只需忽略数据集中 x 值为 None 的行即可。

class MyTransformer(dc.trans.Transformer):
def transform_array(x, y, w, ids):
    kept_rows = x != None
    return x[kept_rows], y[kept_rows], w[kept_rows], ids[kept_rows],

datasets = [d.transform(MyTransformer) for d in datasets]

print(datasets)

train, val, test = datasets

训练模型#

现在我们已经有了数据集，让我们继续并初始化一个 AtomicConvModel 来进行训练。保持输入参数与 AtomicConvFeaturizer 中使用的相同，否则我们将得到错误。 layer_sizes 控制层数和网络中每个密集层的大小。我们选择这些超参数与原文中使用的相同。

acm = AtomicConvModel(n_tasks=1,
                    frag1_num_atoms=f1_num_atoms,
                    frag2_num_atoms=f2_num_atoms,
                    complex_num_atoms=f1_num_atoms+f2_num_atoms,
                    max_num_neighbors=max_num_neighbors,
                    batch_size=12,
                    layer_sizes=[32, 32, 16],
                    learning_rate=0.003,
                    )

losses, val_losses = [], []

%%time
max_epochs = 50

metric = dc.metrics.Metric(dc.metrics.score_function.rms_score)
step_cutoff = len(train)//12
def val_cb(model, step):
if step%step_cutoff!=0:
    return
val_losses.append(model.evaluate(val, metrics=[metric])['rms_score']**2)  # L2 Loss
losses.append(model.evaluate(train, metrics=[metric])['rms_score']**2)  # L2 Loss

acm.fit(train, nb_epoch=max_epochs, max_checkpoints_to_keep=1,
                callbacks=[val_cb])

损失曲线并不是完全平滑的，这并不奇怪，因为我们使用了154个训练和19个验证数据点。增加数据集的大小可能会有所帮助，但也需要更多的计算资源。

f, ax = plt.subplots()
ax.scatter(range(len(losses)), losses, label='train loss')
ax.scatter(range(len(val_losses)), val_losses, label='val loss')
plt.legend(loc='upper right')

ACNN 论文显示，pbbind 核心训练/测试集随机 80/20 拆分的 Pearson \(R^2\) 分数为 0.912 和 0.448。在这里，我们使用 80/10/10 训练/验证/测试拆分，并在训练集（0.943）上实现了类似的性能。我们可以从训练、验证和测试集的表现（以及从论文的结果）看出，ACNN 可以从小型训练数据集学习化学相互作用，但难以一般化。尽管如此，我们只用几行代码就可以训练一个 AtomicConvModel ，并开始预测结合亲和力，这是非常令人惊讶的！

从这里开始，你可以尝试不同的超参数、更具挑战性的拆分和“精炼（refined）”的 PDBbind 数据集，以查看是否可以减少过拟合，并提出一个更严谨的模型。

score = dc.metrics.Metric(dc.metrics.score_function.pearson_r2_score)
for tvt, ds in zip(['train', 'val', 'test'], datasets):
    print(tvt, acm.evaluate(ds, metrics=[score]))

延伸阅读#

我们已经探索了 ACNN 的架构，并使用 PDBbind 数据集训练 ACNN 预测蛋白质-配体结合能。要了解更多信息，请阅读介绍 ACNN 的原始论文：Gomes, Joseph, et al. “Atomic convolutional networks for predicting protein-ligand binding affinity.” arXiv preprint arXiv:1703.10603 (2017)。在预测结合亲和力方面还有许多其他的方法和论文。这里有一些有趣的例子：只用配体或蛋白质进行预测，使用深度学习的分子对接和 AtomNet 。

有条件的生成对抗网络（Conditional Generative Adversarial Network）#

背景知识：

numpy.einsum：推荐观看此视频

生成对抗网络（GAN）是一种生成模型。它由两个部分组成，分别是“生成器”和“判别器”。生成器接受随机值作为输入，并将它们转换为（希望）类似于训练数据的输出。判别器接受一组样本作为输入，并试图将真正的训练样本与生成器创建的样本区分开来。他们都是一起接受训练的。判别器试图越来越好地区分真实数据和虚假数据，而生成器则试图越来越好地欺骗判别器。

有条件的 GAN （CGAN）允许向生成器和判别器输入附加的输出条件。例如，这可能是一个类标签，GAN 试图了解类之间的数据分布如何变化。

在本例中，我们将创建一个由一组分布在椭圆中的二维数据，每个椭圆的位置、形状和方向都是随机的。每个类对应于一个不同的椭圆。让我们随机生成椭圆。对于每一个，我们选择一个随机的中心位置，X 和 Y 的大小和旋转角度。然后我们创建一个变换矩阵将单位圆映射到椭圆。（具体做法，我也没搞清楚，反正下面代码会生成四组shape为(1000,2)的数据集，每一组数据集都是在一个特定椭圆的范围内，四组数据集代表着四个椭圆，这四个数据集保存在points中，也就是说points的shape为(4,1000,2)，这四个椭圆是随机生成的）

import deepchem as dc
import numpy as np
import tensorflow as tf

n_classes = 4
class_centers = np.random.uniform(-4, 4, (n_classes, 2))
class_transforms = []
for i in range(n_classes):
    xscale = np.random.uniform(0.5, 2)
    yscale = np.random.uniform(0.5, 2)
    angle = np.random.uniform(0, np.pi)
    m = [[xscale*np.cos(angle), -yscale*np.sin(angle)],
        [xscale*np.sin(angle), yscale*np.cos(angle)]]
    class_transforms.append(m)
class_transforms = np.array(class_transforms)

该函数从分布中生成随机数据。对于每个点，它随机选择一个类，然后在该类椭圆中随机选择一个位置。

def generate_data(n_points):
    classes = np.random.randint(n_classes, size=n_points)
    r = np.random.random(n_points)
    angle = 2*np.pi*np.random.random(n_points)
    points = (r*np.array([np.cos(angle), np.sin(angle)])).T
    points = ('ijk,ik->ij', class_transforms[classes], points)
    points += class_centers[classes]
    return classes, points

我们从这个分布中画出一些随机点来看看它是什么样子的。点的颜色基于它们的类标签。（从图中可以看出，四个不同的颜色大致构成四个不同的椭圆）

%matplotlib inline
import matplotlib.pyplot as plot
classes, points = generate_data(1000)
plot.scatter(x=points[:,0], y=points[:,1], c=classes)

现在让我们为 CGAN 创建模型。DeepChem 的 GAN 类使这变得非常容易。我们只是创建它的子类并编程一些方法。其中最重要的两个是：

create_generator() 构建一个实现生成器的模型。该模型将一批随机噪声加上任何条件变量（在我们的例子中，每个样本的 one-hot 编码类）作为输入。它的输出是一个合成样本，应该与训练数据相似。

create_discriminator() 构建一个实现判别器的模型。该模型将要评估的样本（可以是真实的训练数据，也可以是生成器创建的合成样本）和条件变量作为输入。它的输出是对应于每个样本的单个数字，这将被解释为样本是真实训练数据的概率。

在本例中，我们使用非常简单的模型。它们只是将输入连接在一起，并通过一些密集层。注意，判别器的最后一层使用了 sigmoid 激活函数。这确保它产生的输出在0到1之间，可以被解释为概率。

我们还需要编码一些定义各种输入形状的方法。我们规定提供给生成器的每个随机噪声样本应该由10个数字组成；每个数据样本包含两个数字（对应一个点的二维X坐标和Y坐标）；条件输入由每个样本的 n_classes 数字组成（one-hot 编码类索引）。

from tensorflow.keras.layers import Concatenate, Dense, Input

class ExampleGAN(dc.models.GAN):

def get_noise_input_shape(self):
    return (10,)

def get_data_input_shapes(self):
    return [(2,)]

def get_conditional_input_shapes(self):
    return [(n_classes,)]

def create_generator(self):
    noise_in = Input(shape=(10,))
    conditional_in = Input(shape=(n_classes,))
    gen_in = Concatenate()([noise_in, conditional_in])
    gen_dense1 = Dense(30, activation=tf.nn.relu)(gen_in)
    gen_dense2 = Dense(30, activation=tf.nn.relu)(gen_dense1)
    generator_points = Dense(2)(gen_dense2)
    return tf.keras.Model(inputs=[noise_in, conditional_in], outputs=[generator_points])

def create_discriminator(self):
    data_in = Input(shape=(2,))
    conditional_in = Input(shape=(n_classes,))
    discrim_in = Concatenate()([data_in, conditional_in])
    discrim_dense1 = Dense(30, activation=tf.nn.relu)(discrim_in)
    discrim_dense2 = Dense(30, activation=tf.nn.relu)(discrim_dense1)
    discrim_prob = Dense(1, activation=tf.sigmoid)(discrim_dense2)
    return tf.keras.Model(inputs=[data_in, conditional_in], outputs=[discrim_prob])

gan = ExampleGAN(learning_rate=1e-4)

来拟合模型。我们通过调用 fit_gan() 来做到这一点。实参是一个生成批训练数据集的迭代器。更具体地说，它需要生成字典，将所有数据输入和条件输入映射到用于它们的值。在我们的例子中，我们可以很容易地创建任意数量的随机数据，因此我们定义了一个生成器，它为每批调用上面定义的 generate_data() 函数。

def iterbatches(batches):
    for i in range(batches):
        classes, points = generate_data(gan.batch_size)
        classes = dc.metrics.to_one_hot(classes, n_classes)
        yield {gan.data_inputs[0]: points, gan.conditional_inputs[0]: classes}

gan.fit_gan(iterbatches(5000))

让训练过的模型生成一些数据，看看它与我们之前绘制的训练分布匹配得有多好。

classes, points = generate_data(1000)
one_hot_classes = dc.metrics.to_one_hot(classes, n_classes)
gen_points = gan.predict_gan_generator(conditional_inputs=[one_hot_classes])
plot.scatter(x=gen_points[:,0], y=gen_points[:,1], c=classes)

完。

在MNIST数据集上训练一个生成对抗网络#

在本教程中，我们将在 MNIST 数据集上训练一个生成对抗网络（GAN）。这是一个手写数字的28x28像素图像的大集合。我们将尝试训练一个网络来生成新的手写数字图像。

首先，让我们导入所有需要的库并加载数据集（它与 Tensorflow 捆绑在一起）。

import deepchem as dc
import tensorflow as tf
from deepchem.models.optimizers import ExponentialDecay
from tensorflow.keras.layers import Conv2D, Conv2DTranspose, Dense, Reshape
import matplotlib.pyplot as plot
import matplotlib.gridspec as gridspec
%matplotlib inline

mnist = tf.keras.datasets.mnist.load_data(path='mnist.npz')
images = mnist[0][0].reshape((-1, 28, 28, 1))/255
dataset = dc.data.NumpyDataset(images)

让我们来看看其中的一些图片，了解一下它们的样子。

def plot_digits(im):
plot.figure(figsize=(3, 3))
grid = gridspec.GridSpec(4, 4, wspace=0.05, hspace=0.05)
for i, g in enumerate(grid):
    ax = plot.subplot(g)
    ax.set_xticks([])
    ax.set_yticks([])
    ax.imshow(im[i,:,:,0], cmap='gray')

plot_digits(images)

现在我们可以创建 GAN 了。和上一篇教程一样，它由两部分组成:

生成器会把随机噪声作为输入，并试图产生与训练数据相似的输出。
判别器接受一组样本作为输入（可能是训练数据，可能是生成器创建的），并试图确定哪些是哪些。

这次我们将使用一种不同风格的 GAN，称为 Wasserstein GAN （或简称WGAN）。在许多情况下，它们被发现比传统的 GANs 产生更好的结果。两者之间的主要区别在于判别器（在这种情况下通常称为“批评家”）。它试图学习如何测量训练数据集分布和生成数据集分布之间的距离，而不是输出样本作为真实训练数据的概率。然后，该测量方法可以直接用来作为训练生成器的损失函数。

我们使用一个非常简单的模型。该生成器使用一个密集层将噪声输入转换为一个具有8通道的7x7图像。接下来是两个卷积层，首先升采样到14x14，最后升采样到28x28。

判别器反向执行大致相同的操作。两个卷积层降采样图像，首先到14x14，然后到7x7。最后一个密集层产生一个数字作为输出。在上一篇教程中，我们使用了 sigmoid 激活来生成一个0到1之间的数字，它可以被解释为概率。由于这是一个 WGAN ，因此我们使用 softplus 激活函数。它产生一个无界正数，可以解释为一个距离。

class DigitGAN(dc.models.WGAN):

def get_noise_input_shape(self):
    return (10,)

def get_data_input_shapes(self):
    return [(28, 28, 1)]

def create_generator(self):
    return tf.keras.Sequential([
        Dense(7*7*8, activation=tf.nn.relu),
        Reshape((7, 7, 8)),
        Conv2DTranspose(filters=16, kernel_size=5, strides=2, activation=tf.nn.relu, padding='same'),
        Conv2DTranspose(filters=1, kernel_size=5, strides=2, activation=tf.sigmoid, padding='same')
    ])

def create_discriminator(self):
    return tf.keras.Sequential([
        Conv2D(filters=32, kernel_size=5, strides=2, activation=tf.nn.leaky_relu, padding='same'),
        Conv2D(filters=64, kernel_size=5, strides=2, activation=tf.nn.leaky_relu, padding='same'),
        Dense(1, activation=tf.math.softplus)
    ])

gan = DigitGAN(learning_rate=ExponentialDecay(0.001, 0.9, 5000))

现在训练它。和上一篇教程一样，我们编写了一个数据生成器来生成数据。这次的数据来自一个数据集，我们将循环100次。

还有一个区别值得注意。在训练传统 GAN 时，重要的是在整个训练过程中保持生成器和判别器的平衡。如果其中一个走得太远，另一个就很难学习了。

WGANs 没有这个问题。事实上，判别器变得越好，它提供的信号就越清晰，生成器就越容易学习它。因此，我们指定 generator_steps=0.2 ，以便每训练判别器5步，只需要训练生成器1步。这往往产生更快的训练效果和更好的结果。

def iterbatches(epochs):
    for i in range(epochs):
        for batch in dataset.iterbatches(batch_size=gan.batch_size):
        yield {gan.data_inputs[0]: batch[0]}
    gan.fit_gan(iterbatches(100), generator_steps=0.2, checkpoint_interval=5000)

让我们生成一些数据，看看结果如何。

plot_digits(gan.predict_gan_generator(batch_size=16))

还不错。许多生成的图像看起来很像手写的数字。当然，经过较长时间训练的较大模型可以做得更好。

LitMatter DeepChem#

Jupyter Notebook 需要下载整个 Github代码然后再从 litmatter 文件夹中打开 LitDeepChem.ipynb

git clone https://github.com/ncfrey/litmatter.git

Jupyter Notebook 中文翻译版下载需要下载整个 Github 代码然后再从 litmatter 文件夹中打开 LitDeepChem.ipynb

git clone https://github.com/abdusemiabduweli/AIDD-Tutorial-Files.git

这本笔记本展示了如何使用 LitMatter 模板在 MoleculeNet 数据集上加速 DeepChem 模型训练。
在本例中，我们在 Tox21 数据集上训练一个简单的 DeepChem TorchModel 。
这里展示的训练工作流可以通过更改一个关键参数扩展到数百个 GPU！

%load_ext autoreload
%autoreload 2

import torch

import deepchem as dc

import pytorch_lightning as pl
from pytorch_lightning.callbacks import ModelCheckpoint
from pytorch_lightning import (LightningDataModule, LightningModule, Trainer,
                            seed_everything)

加载一个 LitMolNet 数据集#

任何来自 deepchem.molnet 的 MolNet 数据集可与 LitMatter 配合使用。具体的 MolNet 数据集和任何预处理步骤都可以在 data.LitMolNet 中定义。

from lit_data.molnet_data import LitMolNet

dm = LitMolNet(loader=dc.molnet.load_tox21, batch_size=16)
dm.prepare_data()
dm.setup()

实例化一个 LitDeepChem 模型#

任何 deepchem.models.torch_models.TorchModel 可以与 LitMatter 一起使用。在这里，我们将在 PyTorch 中编写我们自己的自定义基本模型，并创建一个 TorchModel 。

from lit_models.deepchem_models import LitDeepChem

base_model = torch.nn.Sequential(
torch.nn.Linear(1024, 256),
    torch.nn.ReLU(),
    torch.nn.Linear(256, 12),
)

torch_model = dc.models.TorchModel(base_model, loss=torch.nn.MSELoss())

model = LitDeepChem(torch_model, lr=1e-2)

训练模型#

在用多个 GPU 和多节点训练时，只需根据需要更改 Trainer 标志。

trainer = Trainer(gpus=-1,  # use all available GPUs on each node
#                   num_nodes=1,  # change to number of available nodes
#                  accelerator='ddp',
                max_epochs=5,
                )

trainer.fit(model, datamodule=dm)

就这样！通过改变 num_nodes 参数，训练可以分布在所有可用的 GPU 上。有关 HPC 集群上较长的训练作业，请参阅提供的示例批处理脚本。

使用 Hyperopt 高级模型训练#

在高级模型训练教程中，我们已经了解了在 deepchem 包中使用 GridHyperparamOpt 进行超参数优化。在本教程中，我们将研究另一个称为 Hyperopt 的超参数调优库。

准备#

要运行本教程需要安装 Hyperopt 库。

pip install hyperopt

通过hyperopt进行超参数优化#

让我们从加载 HIV 数据集开始。它根据是否抑制艾滋病毒复制对超过4万个分子进行了分类。

import deepchem as dc
tasks, datasets, transformers = dc.molnet.load_hiv(featurizer='ECFP', split='scaffold')
train_dataset, valid_dataset, test_dataset = datasets

现在，让我们导入 hyperopt 库，我们将使用它来提供最佳参数

from hyperopt import hp, fmin, tpe, Trials

然后，我们必须声明一个字典，其中包含所有超形参及其将调优的范围。这本字典将作为 hyperopt 的搜索空间。

在字典中声明范围的一些基本方法是：

hp.choice(‘label’,[choices]) : this is used to specify a list of choices
hp.uniform(‘label’ ,low= low_value ,high= high_value ) : this is used to specify a uniform distibution

在低值和高值之间。它们之间的值可以是任何实数，不一定是整数。

在这里，我们将使用多任务分类器对 HIV 数据集进行分类，因此适当的搜索空间如下所示。

search_space = {
    'layer_sizes': hp.choice('layer_sizes',[[500], [1000], [2000],[1000,1000]]),
    'dropouts': hp.uniform('dropout',low=0.2, high=0.5),
    'learning_rate': hp.uniform('learning_rate',high=0.001, low=0.0001)
}

然后，我们应该声明一个由 hyperopt 最小化的函数。所以，这里我们应该使用这个函数来最小化我们的多任务分类器模型。此外，我们使用 validation callback 在每1000步验证分类器，然后将最佳分数作为返回值传递。这里使用的指标是 ` roc_auc_score ` ，需要最大化它。使一个非负值最大化相当于使其负数最小化，因此我们将返回验证得分的负数。

import tempfile
#tempfile is used to save the best checkpoint later in the program.

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

def fm(args):
save_dir = tempfile.mkdtemp()
model = dc.models.MultitaskClassifier(n_tasks=len(tasks),n_features=1024,layer_sizes=args['layer_sizes'],dropouts=args['dropouts'],learning_rate=args['learning_rate'])
#validation callback that saves the best checkpoint, i.e the one with the maximum score.
validation=dc.models.ValidationCallback(valid_dataset, 1000, [metric],save_dir=save_dir,transformers=transformers,save_on_minimum=False)

model.fit(train_dataset, nb_epoch=25,callbacks=validation)

#restoring the best checkpoint and passing the negative of its validation score to be minimized.
model.restore(model_dir=save_dir)
valid_score = model.evaluate(valid_dataset, [metric], transformers)

return -1*valid_score['roc_auc_score']

在这里，我们调用 hyperopt 的 fmin 函数，在这里我们传递要最小化的函数、要遵循的算法、最大 eval 数和一个 trials 对象。Trials 对象用于保存所有超参数、损失和其他信息，这意味着你可以在运行优化后访问它们。此外，Trials 可以帮助你保存重要信息，以便稍后加载，然后恢复优化过程。

此外，该算法有三种选择，无需额外配置即可使用。他们是:-

Random Search - rand.suggest
TPE (Tree Parzen Estimators) - tpe.suggest
Adaptive TPE - atpe.suggest

下面的代码用于打印 hyperopt 找到的最佳超参数。

print("Best: {}".format(best))

这里发现的超参数不一定是最好的，但可以大致了解哪些参数是有效的。为了得到更准确的结果，必须增加验证周期的数量和模型拟合的周期。但是这样做可能会增加寻找最佳超参数的时间。

高斯过程（Gaussian Processes）简介#

在化学信息学和机器学习的世界中，模型通常是树（随机森林、XGBoost 等）或人工神经网络（深度神经网络、图卷积网络等）。这些模型被称为“频率主义者（Frequentist）”模型。然而，还有另一类被称为贝叶斯模型。今天我们将实验一个在 scikit-learn 中实现的贝叶斯模型，称为高斯过程（gaussian processes，GP）。要深入了解GP，有一个很好的教程论文介绍了GP如何用于回归。还有一篇学术论文将GP应用于一个现实世界的问题。

作为一个简短的介绍，GP 允许我们在 n 维空间上使用无限个高斯函数建立统计模型，其中 n 是特征的数量。然而，我们选择这些函数是基于它们对传递给它的数据的拟合程度。我们最终得到了一个由不固定的高斯函数的集合建立的统计模型。结果是，对于我们训练模型的点，在我们的集合中，方差应该非常低。对于接近训练集合点的测试集合点，方差应该较高，但仍然较低，因为函数集合为了在其邻近区域内预测良好而选取的。然而，对于远离训练集合点的点，我们没有选择可以拟合它们的高斯函数的集合，因此我们预计在函数集合中方差会很高。通过这种方式，我们最终得到了一个统计模型，它允许不确定性的自然产生。

高斯过程#

如前所述，GP 已经在 scikit-learn 中实现，因此我们将使用 DeepChem 的 scikit-learn 包装器。SklearnModel 是 DeepChem Model 类的子类。它充当 sklearn.base.BaseEstimator 的包装器。

这里我们导入 deepchem 和从 sklearn 中 GP 回归模型。

import deepchem as dc
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel

import numpy as np
import matplotlib.pyplot as plt

加载数据#

接下来，我们需要一个反映回归问题的数据集。对于本教程，我们将使用来自 MoleculeNet 的 BACE 数据集。

tasks, datasets, transformers = dc.molnet.load_bace_regression(featurizer='ecfp', splitter='random')
train_dataset, valid_dataset, test_dataset = datasets

我总是喜欢仔细查看代码中的对象存储了什么。我们看到任务是我们试图预测的一系列任务。转换器是 NormalizationTransformer，它规范化数据集的输出（y值）。

print(f'The tasks are: {tasks}')
print(f'The transformers are: {transformers}')
print(f'The transformer normalizes the outputs (y values): {transformers[0].transform_y}')

在这里，我们可以看到数据已经被分成了一个训练集、一个验证集和一个测试集。我们将在训练集中训练模型，并在测试集中测试模型的准确性。如果要进行任何超参数调优，我们将使用验证集。拆分是大约80/10/10次训练/验证/测试。

print(train_dataset)
print(valid_dataset)
print(test_dataset)

使用 SklearnModel#

这里，我们首先使用从 sklearn 导入的 GaussianProcessRegressor 创建模型。然后我们把它包装在 DeepChem 的 SklearnModel 中。要了解关于模型的更多信息，可以阅读 sklearn API 或在代码块中运行 help(GaussianProcessRegressor) 。

如你所见，我为参数选择的值似乎非常具体。这是因为我需要事先做一些超参数调优，以得到不会过度拟合训练集的模型。你可以在本教程末尾的附录中了解更多关于我如何调优模型的信息。

output_variance = 7.908735015054668
length_scale = 6.452349252677817
noise_level = 0.10475507755839343
kernel = output_variance**2 * RBF(length_scale=length_scale, length_scale_bounds='fixed') + WhiteKernel(noise_level=noise_level, noise_level_bounds='fixed')
alpha = 4.989499481123432e-09

sklearn_gpr = GaussianProcessRegressor(kernel=kernel, alpha=alpha)
model = dc.models.SklearnModel(sklearn_gpr)

然后我们将模型拟合到数据，看看它在训练集和测试集上的表现如何。

model.fit(train_dataset)
metric1 = dc.metrics.Metric(dc.metrics.mean_squared_error)
metric2 = dc.metrics.Metric(dc.metrics.r2_score)
print(f'Training set score: {model.evaluate(train_dataset, [metric1, metric2])}')
print(f'Test set score: {model.evaluate(test_dataset, [metric1, metric2])}')

分析结果#

我们还可以将预测值与实测值的匹配程度可视化。首先，我们需要一个函数，它允许我们获得预测值的平均值和该值的标准差。这是通过从每组输入 X 中抽样100个预测，并计算平均值和标准差来实现的。

def predict_with_error(dc_model, X, y_transformer):
    samples = model.model.sample_y(X, 100)
    means = y_transformer.untransform(np.mean(samples, axis=1))
    stds = y_transformer.y_stds[0] * np.std(samples, axis=1)

    return means, stds

对于我们的训练集，我们看到测量值(x轴)和预测值(y轴)之间有很好的相关性。注意，我们使用前面的转换器来取消对预测值的转换。

y_meas_train = transformers[0].untransform(train_dataset.y)
y_pred_train, y_pred_train_stds = predict_with_error(model, train_dataset.X, transformers[0])

plt.xlim([2.5, 10.5])
plt.ylim([2.5, 10.5])
plt.scatter(y_meas_train, y_pred_train)

现在，我们对我们的测试集执行同样的操作。我们看到了相当好的相关性！然而，它肯定没有那么紧。这反映在上面计算的 R2 分数之间的差异上。

y_meas_test = transformers[0].untransform(test_dataset.y)
y_pred_test, y_pred_test_stds = predict_with_error(model, test_dataset.X, transformers[0])

plt.xlim([2.5, 10.5])
plt.ylim([2.5, 10.5])
plt.scatter(y_meas_test, y_pred_test)

我们还可以编写一个函数来计算有多少预测值落在预测误差范围内。这是通过计算有多少个样本的真实误差小于之前计算的标准偏差来实现的。一个标准差是68%的置信区间。

def percent_within_std(y_meas, y_pred, y_std):
    assert len(y_meas) == len(y_pred) and len(y_meas) == len(y_std), 'length of y_meas and y_pred must be the same'

    count_within_error = 0
    for i in range(len(y_meas)):
        if abs(y_meas[i][0]-y_pred[i]) < y_std[i]:
            count_within_error += 1

    return count_within_error/len(y_meas)

对于训练集，90%的样本在一个标准偏差内。相比之下，只有约70%的样本在测试集的标准偏差内。一个标准差是68%的置信区间所以我们看到对于训练集，不确定度很接近。然而，该模型过度预测了训练集的不确定性。

print(percent_within_std(y_meas_train, y_pred_train, y_pred_train_stds))
print(percent_within_std(y_meas_test, y_pred_test, y_pred_test_stds))

我们还可以看看测试集预测的标准偏差的分布。我们在预测误差中看到一个非常粗略的高斯分布。

plt.hist(y_pred_test_stds)
plt.show()

现在，我们的教程到此结束。我们计划很快继续深入研究不确定性估计，特别是校准的不确定性估计。我们到时候见！

附录：Hyperparameter 优化#

由于超参数优化超出了本教程的范围，所以我将不解释如何使用Optuna调优超参数。但是为了完整起见，仍然包含了代码。

%pip install optuna
import optuna


def get_model(trial):
    output_variance = trial.suggest_float('output_variance', 0.1, 10, log=True)
    length_scale = trial.suggest_float('length_scale', 1e-5, 1e5, log=True)
    noise_level = trial.suggest_float('noise_level', 1e-5, 1e5, log=True)

    params = {
        'kernel': output_variance**2 * RBF(length_scale=length_scale, length_scale_bounds='fixed') + WhiteKernel(noise_level=noise_level, noise_level_bounds='fixed'),
        'alpha': trial.suggest_float('alpha', 1e-12, 1e-5, log=True),
    }

    sklearn_gpr = GaussianProcessRegressor(**params)
    return dc.models.SklearnModel(sklearn_gpr)

def objective(trial):
    model = get_model(trial)
    model.fit(train_dataset)

    metric = dc.metrics.Metric(dc.metrics.mean_squared_error)
    return model.evaluate(valid_dataset, [metric])['mean_squared_error']

study = optuna.create_study(direction='minimize')
study.optimize(objective, n_trials=100)

print(study.best_params)

与PytorchLightning结合#

在本教程中，我们将介绍如何在 pytorch-lightning 框架中设置 deepchem 模型。Lightning 是一个 pytorch 框架，它简化了使用 pytorch 模型的实验过程。pytorch lightning 提供的以下几个关键功能是 deepchem 用户可以发现有用的:

多 gpu 训练功能：pytorch-lightning 提供简单的多 gpu、多节点训练。它还简化了跨不同集群基础设施（如 AWS、基于 slurm 的集群）启动多gpu、多节点作业的过程。
减少 pytorch 的样板代码：lightning 负责处理诸如 optimizer.zero_grad(), model.train(), model.eval() 之类的细节。Lightning 还提供了实验日志功能，例如，无论在 CPU、GPU、多节点上进行训练，用户都可以在训练器内部使用 self.log 方法，它将适当地记录指标。
可以加速训练的工具：半精确训练、梯度检查点、代码分析。

准备#

本笔记本假设你已经安装了deepchem，如果你没有，请执行 deepchem 安装页面的说明：https://deepchem.readthedocs.io/en/latest/get_started/installation.html。
安装 pytorchlightning 请参考 lightning 的主页：https://www.pytorchlightning.ai/

导入相关的包。

import deepchem as dc
from deepchem.models import GCNModel
import pytorch_lightning as pl
import torch
from torch.nn import functional as F
from torch import nn
import pytorch_lightning as pl
from pytorch_lightning.core.lightning import LightningModule
from torch.optim import Adam
import numpy as np
import torch

Deepchem 例子#

下面我们展示图卷积网络（GCN）的一个例子。请注意，这是一个简单的示例，它使用 GCNModel 从输入序列预测标签。在这个例子中，我们没有展示 deepchem 的完整功能，因为我们想要重组 deepchem 代码，并对其进行调整，以便能够轻松插入 pytorch-lightning。这个例子的灵感来自 GCNModel 文档。

准备数据集：为了训练我们的deepchem模型，我们需要一个可以用来训练模型的数据集。下面我们为本教程准备了一个示例数据集。下面我们也直接使用特征器编码数据集编码。

smiles = ["C1CCC1", "CCC"]
labels = [0., 1.]
featurizer = dc.feat.MolGraphConvFeaturizer()
X = featurizer.featurize(smiles)
dataset = dc.data.NumpyDataset(X=X, y=labels)

设置模型：现在我们初始化我们将在训练中使用的图卷积网络模型。

model = GCNModel(
    mode='classification',
    n_tasks=1,
    batch_size=2,
    learning_rate=0.001
)

训练模型：在我们的训练数据集上拟合模型，还指定要运行的epoch的数量。

loss = model.fit(dataset, nb_epoch=5)
print(loss)

Pytorch-Lightning + Deepchem 示例#

现在我们来看一个 GCN 模型适用于 Pytorch-Lightning 的例子。使用 Pytorch-Lightning 有两个重要的组成部分：

LightningDataModule ：该模块定义如何准备数据并将其输入到模型中，以便模型可以使用它进行训练。该模块定义了训练数据加载器函数，训练器直接使用该函数为 LightningModule 生成数据。要了解有关 LightningDataModule 的更多信息，请参阅 datamodules 文档。
LightningModule ：这个模块为我们的模型定义了训练和验证步骤。我们可以使用这个模块根据超参数初始化我们的模型。我们可以直接使用许多样板函数来跟踪我们的实验，例如，我们可以使用 self.save_hyperparameters() 方法来保存我们用于训练的所有超参数。有关如何使用该模块的详细信息，请参阅 lightningmodules 文档。

设置 torch 数据集：请注意，这里我们需要创建一个自定义的 SmilesDataset，以便我们可以轻松地与 deepchem 特征器交互。为这个交互,我们需要定义一个整理方法,这样我们就可以创建批次数据集。

# prepare LightningDataModule
class SmilesDataset(torch.utils.data.Dataset):
    def __init__(self, smiles, labels):
        assert len(smiles) == len(labels)
        featurizer = dc.feat.MolGraphConvFeaturizer()
        X = featurizer.featurize(smiles)
        self._samples = dc.data.NumpyDataset(X=X, y=labels)

    def __len__(self):
        return len(self._samples)

    def __getitem__(self, index):
        return (
            self._samples.X[index],
            self._samples.y[index],
            self._samples.w[index],
        )


class SmilesDatasetBatch:
    def __init__(self, batch):
        X = [np.array([b[0] for b in batch])]
        y = [np.array([b[1] for b in batch])]
        w = [np.array([b[2] for b in batch])]
        self.batch_list = [X, y, w]


def collate_smiles_dataset_wrapper(batch):
    return SmilesDatasetBatch(batch)

创建 lightning 模块：在本部分中，我们创建 GCN 特定的 lightning 模块。该类指定训练步骤的逻辑流。我们还为训练流创建所需的模型、优化器和损失。

# prepare the LightningModule
class GCNModule(pl.LightningModule):
    def __init__(self, mode, n_tasks, learning_rate):
        super().__init__()
        self.save_hyperparameters(
            "mode",
            "n_tasks",
            "learning_rate",
        )
        self.gcn_model = GCNModel(
            mode=self.hparams.mode,
            n_tasks=self.hparams.n_tasks,
            learning_rate=self.hparams.learning_rate,
        )
        self.pt_model = self.gcn_model.model
        self.loss = self.gcn_model._loss_fn

    def configure_optimizers(self):
        return self.gcn_model.optimizer._create_pytorch_optimizer(
            self.pt_model.parameters(),
        )

    def training_step(self, batch, batch_idx):
        batch = batch.batch_list
        inputs, labels, weights = self.gcn_model._prepare_batch(batch)
        outputs = self.pt_model(inputs)

        if isinstance(outputs, torch.Tensor):
            outputs = [outputs]

        if self.gcn_model._loss_outputs is not None:
            outputs = [outputs[i] for i in self.gcn_model._loss_outputs]

        loss_outputs = self.loss(outputs, labels, weights)

        self.log(
            "train_loss",
            loss_outputs,
            on_epoch=True,
            sync_dist=True,
            reduce_fx="mean",
            prog_bar=True,
        )

        return loss_outputs

创建相关对象

# create module objects
smiles_datasetmodule = SmilesDatasetModule(
    train_smiles=["C1CCC1", "CCC", "C1CCC1", "CCC", "C1CCC1", "CCC", "C1CCC1", "CCC", "C1CCC1", "CCC"],
    train_labels=[0., 1., 0., 1., 0., 1., 0., 1., 0., 1.],
    batch_size=2,
)

gcnmodule = GCNModule(
    mode="classification",
    n_tasks=1,
    learning_rate=1e-3,
)

Lightning 训练器#

Trainer 是构建在 LightningDataModule 和 LightningModule 之上的包装器。当构建 lightning 训练器时，你还可以指定 epoch 的数量，运行的最大步数，gpu 的数量，用于训练器的节点数量。Lightning trainer 充当分布式训练设置的包装器，这样你就能够简单地构建模型以本地运行。

trainer = pl.Trainer(
    max_epochs=5,
)

调用 fit 函数运行模型训练

# train
trainer.fit(
    model=gcnmodule,
    datamodule=smiles_datasetmodule,
)

分子无监督嵌入学习#

在本教程中，我们将使用 SeqToSeq 模型来生成用于分类分子的指纹。这是基于以下论文，尽管一些实现细节是不同的：Xu et al., “Seq2seq Fingerprint: An Unsupervised Deep Molecular Embedding for Drug Discovery” (https://doi.org/10.1145/3107411.3107424).

使用SeqToSeq学习嵌入#

许多类型的模型要求它们的输入具有固定的形状。由于分子所包含的原子和化学键的数量差异很大，因此很难将这些模型应用到分子上。我们需要一种方法为每个分子生成固定长度的“指纹”。为此设计了各种方法，例如我们在前面的教程中使用的扩展连接指纹(Extended-Connectivity prints, ECFPs)。但在这个例子中，我们将让 SeqToSeq 模型学习自己创建指纹的方法，而不是手工设计指纹。

SeqToSeq 模型执行序列到序列的翻译。例如，它们经常被用来将文本从一种语言翻译成另一种语言。它由“编码器”和“解码器”两部分组成。编码器是一个循环层的堆栈。输入序列被输入到它中，一次一个token，它生成一个固定长度的向量，称为“嵌入向量”。解码器是执行反向操作的另一个循环层堆栈:它接受嵌入向量作为输入，并生成输出序列。通过对适当选择的输入/输出对进行训练，你可以创建一个执行多种转换的模型。

在本例中，我们将使用描述分子的SMILES字符串作为输入序列。我们将把模型训练成一个自动编码器，这样它就会尝试使输出序列与输入序列相同。为此，编码器必须创建包含原始序列中所有信息的嵌入向量。这正是我们在指纹中想要的，所以也许这些嵌入向量将会在其他模型中作为一种表示分子的方式有用！

让我们从加载数据开始。我们将使用MUV数据集。它在训练集中包含74501个分子，在验证集中包含9313个分子，因此它给我们提供了大量的SMILES字符串来处理。

import deepchem as dc
tasks, datasets, transformers = dc.molnet.load_muv(split='stratified')
train_dataset, valid_dataset, test_dataset = datasets
train_smiles = train_dataset.ids
valid_smiles = valid_dataset.ids

我们需要为 SeqToSeq 模型定义“字母表”，即可以出现在序列中的所有表示的列表。（输入和输出序列也可能具有不同的字母，但由于我们将其训练为自动编码器，所以在这种情况下它们是相同的。）列出在任何训练序列中出现的每个字符。

tokens = set()
for s in train_smiles:
tokens = tokens.union(set(c for c in s))
tokens = sorted(list(tokens))

创建模型并定义要使用的优化方法。在这种情况下，如果我们逐渐降低学习速度，学习效果会更好。我们使用 ExponentialDecay 在每次迭代后将学习率乘以0.9。

from deepchem.models.optimizers import Adam, ExponentialDecay
max_length = max(len(s) for s in train_smiles)
batch_size = 100
batches_per_epoch = len(train_smiles)/batch_size
model = dc.models.SeqToSeq(tokens,
                        tokens,
                        max_length,
                        encoder_layers=2,
                        decoder_layers=2,
                        embedding_dimension=256,
                        model_dir='fingerprint',
                        batch_size=batch_size,
                        learning_rate=ExponentialDecay(0.001, 0.9, batches_per_epoch))

让我们来训练它！ fit_sequences() 的输入是一个生成输入/输出对的生成器。在一个好的 GPU 上，这应该需要几个小时或更少的时间。

def generate_sequences(epochs):
    for i in range(epochs):
        for s in train_smiles:
            yield (s, s)

model.fit_sequences(generate_sequences(40))

让我们看看它作为自动编码器的效果如何。我们将运行验证集中的前500个分子，看看其中有多少被精确复制。

predicted = model.predict_from_sequences(valid_smiles[:500])
    count = 0
for s,p in zip(valid_smiles[:500], predicted):
    if ''.join(p) == s:
        count += 1
print('reproduced', count, 'of 500 validation SMILES strings')

现在我们试着用编码器来生成分子指纹。我们计算训练和验证数据集中所有分子的嵌入向量，并创建新的数据集，这些数据集的特征向量。数据量非常小，我们可以将所有数据都存储在内存中。

import numpy as np
train_embeddings = model.predict_embeddings(train_smiles)
train_embeddings_dataset = dc.data.NumpyDataset(train_embeddings,
                                                train_dataset.y,
                                                train_dataset.w.astype(np.float32),
                                                train_dataset.ids)

valid_embeddings = model.predict_embeddings(valid_smiles)
valid_embeddings_dataset = dc.data.NumpyDataset(valid_embeddings,
                                                valid_dataset.y,
                                                valid_dataset.w.astype(np.float32),
                                                valid_dataset.ids)

为了分类，我们将使用一个具有一个隐藏层的简单全连接网络。

classifier = dc.models.MultitaskClassifier(n_tasks=len(tasks),
                                                    n_features=256,
                                                    layer_sizes=[512])
classifier.fit(train_embeddings_dataset, nb_epoch=10)

看看它的效果如何。计算训练和验证数据集的 ROC AUC。

metric = dc.metrics.Metric(dc.metrics.roc_auc_score, np.mean, mode="classification")
train_score = classifier.evaluate(train_embeddings_dataset, [metric], transformers)
valid_score = classifier.evaluate(valid_embeddings_dataset, [metric], transformers)
print('Training set ROC AUC:', train_score)
print('Validation set ROC AUC:', valid_score)

合成可行性评价#

[1] https://pubs.acs.org/doi/abs/10.1021/acs.jcim.7b00622

合成可行性#

当运行大规模枚举时，合成可行性是一个问题。通常，被列举的分子很难制造，因此不值得检验，即使它们的其他化学性质在计算机模拟中很好。本教程将介绍如何训练 ScScore 模型 [1]。

该模型的思想是在一对分子上进行训练，其中一个分子比另一个分子“更复杂”。然后神经网络就可以做出分数，试图保持分子的这种成对排序。最终的结果是一个可以给出分子相对复杂性的模型。

这篇论文对reaxys中的每一个反应进行了训练，宣称产物比反应物更复杂。由于这个训练集非常昂贵，我们将对任意分子进行训练，如果一个分子的SMILES字符串较长，则该分子声明为更复杂的分子。在现实世界中，你可以使用任何对项目有意义的复杂性度量。

在本教程中，我们将使用 Tox21 数据集来训练我们简单的合成可行性模型。

创建数据集#

让我们从加载一些分子开始。我们加载 Tox21，指定 splitter=None，这样所有内容都将作为单个数据集返回。

import deepchem as dc
tasks, datasets, transformers = dc.molnet.load_tox21(featurizer='Raw', splitter=None)
molecules = datasets[0].X

因为 ScScore 是根据相对复杂度进行训练的，所以我们希望数据集中的 X 张量具有3个维度 (sample_id, molecule_id, features) 。 molecule_id 维度的大小为2，因为一个样本是一对分子。如果第一个分子比第二个分子更复杂，则标记为1。下面我们介绍的函数 create_dataset 从给定的列表中随机提取微笑字符串对，并根据这种复杂性度量对它们进行排序。

在现实世界中，你可以使用购买成本或所需反应步骤数作为复杂度评分。

from rdkit import Chem
import random
from deepchem.feat import CircularFingerprint
import numpy as np


def create_dataset(fingerprints, smiles_lens, ds_size=100000):
    """
    m1: list of np.Array
        fingerprints for molecules
    m2: list of int
        length of a molecules SMILES string

    returns:
        dc.data.Dataset for input into ScScore Model

    Dataset.X
        shape is (sample_id, molecule_id, features)
    Dataset.y
        shape is (sample_id,)
        values is 1 if the 0th index molecule is more complex
                0 if the 1st index molecule is more complex
    """
    X, y = [], []
    all_data = list(zip(fingerprints, smiles_lens))
    while len(y) < ds_size:
        i1 = random.randrange(0, len(smiles_lens))
        i2 = random.randrange(0, len(smiles_lens))
        m1 = all_data[i1]
        m2 = all_data[i2]
        if m1[1] == m2[1]:
            continue
        if m1[1] > m2[1]:
            y.append(1.0)
        else:
            y.append(0.0)
        X.append([m1[0], m2[0]])
    return dc.data.NumpyDataset(np.array(X), np.expand_dims(np.array(y), axis=1))

复杂性排序器就位后，我们现在可以构建数据集了。让我们从随机地将分子列表分成训练集和测试集开始。

molecule_ds = dc.data.NumpyDataset(np.array(molecules))
splitter = dc.splits.RandomSplitter()
train_mols, test_mols = splitter.train_test_split(molecule_ds)

我们将用具有手性的 ECFP 指纹（与源论文匹配）对所有分子进行特征化，然后使用上面定义的函数构建我们的成对数据集。

n_features = 1024
featurizer = dc.feat.CircularFingerprint(size=n_features, radius=2, chiral=True)
train_features = featurizer.featurize(train_mols.X)
train_smiles_len = [len(Chem.MolToSmiles(x)) for x in train_mols.X]
train_dataset = create_dataset(train_features, train_smiles_len)

现在我们已经创建了数据集，让我们在这个数据集上训练一个 ScScoreModel。

model = dc.models.ScScoreModel(n_features=n_features)
model.fit(train_dataset, nb_epoch=20)

模型性能#

让我们来评估一下这个模型在我们的抵制分子上的表现。SaScores 应该从从未见过的分子追踪 SMILES 字符串的长度。

import matplotlib.pyplot as plt
%matplotlib inline

mol_scores = model.predict_mols(test_mols.X)
smiles_lengths = [len(Chem.MolToSmiles(x)) for x in test_mols.X]

现在使用 matplotlib 绘制分子 smiles 字符串的长度与 SaScore 的比较。

plt.figure(figsize=(20,16))
plt.scatter(smiles_lengths, mol_scores)
plt.xlim(0,80)
plt.xlabel("SMILES length")
plt.ylabel("ScScore")
plt.show()

正如我们所看到的，模型通常跟踪 SMILES 的长度。它在8到30个字符之间有很好的丰富，并且大大小小的 SMILES 字符串都非常准确。

现在你可以用比 SMILES 长度更有意义的度量来训练你自己的模型了！

参考文献：#

基于图卷积 QSAR 模型计算原子对分子的贡献#

在前面的教程中，我们介绍了模型可解释性的概念：了解模型为什么会产生它所产生的结果。在本教程中，我们将学习原子对分子的贡献，这是解释作用与分子的模型的有用工具。

想法很简单：从分子中删除一个原子，看看模型的预测如何变化。一个原子的“原子贡献”被定义为整个分子和去除此原子后剩下的碎片之间的活性差。这是一种衡量原子对预测结果影响程度的方法。

贡献在文献中也称为“属性”、“着色”（”attributions”, “coloration”）等。这是一种模型解释方法[1]，类似于 QSAR 域中的 Similarity maps[2]，或其他领域(图像分类等)中的遮挡方法（occlusion methods）。目前的实现在[4]中使用。

https://github.com/deepchem/deepchem/blob/master/examples/tutorials/assets/atomic_contributions_tutorial_data/index.png?raw=1

血脑屏障通透性的 QSAR 分类模型#

血脑屏障的通透性是化合物进入中枢神经系统的能力。在这里，我们使用了一个相对较小的化合物数据集，这些化合物是在没有任何载体的情况下通过扩散运输的。属性定义为log10(大脑浓度/血液浓度)。值为正值(和0)的化合物被标记为活性化合物，其他化合物被标记为非活性化合物。在建模之后，我们将确定有利于扩散和不利于扩散的原子。

首先让我们创建数据集。分子存储在一个 SDF 文件中。

import os
import pandas as pd
import deepchem as dc
import numpy as np
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem import Draw, PyMol, rdFMCS
from rdkit.Chem.Draw import IPythonConsole
from rdkit import rdBase
from deepchem import metrics
from IPython.display import Image, display
from rdkit.Chem.Draw import SimilarityMaps
import tensorflow as tf

current_dir = os.path.dirname(os.path.realpath('__file__'))
dc.utils.download_url(
    'https://raw.githubusercontent.com/deepchem/deepchem/master/examples/tutorials/assets/atomic_contributions_tutorial_data/logBB.sdf',
    current_dir,
    'logBB.sdf'
)
DATASET_FILE =os.path.join(current_dir, 'logBB.sdf')
# Create RDKit mol objects, since we will need them later.
mols = [m for m in Chem.SDMolSupplier(DATASET_FILE) if m is not None ]
loader = dc.data.SDFLoader(tasks=["logBB_class"],
                        featurizer=dc.feat.ConvMolFeaturizer(),
                        sanitize=True)
dataset = loader.create_dataset(DATASET_FILE, shard_size=2000)

现在，让我们构建和训练一个GraphConvModel。

m = dc.models.GraphConvModel(1, mode="classification", batch_normalize=False, batch_size=100)
m.fit(dataset, nb_epoch=10)

让我们加载一个测试集，看看效果如何。

current_dir = os.path.dirname(os.path.realpath('__file__'))
dc.utils.download_url(
    'https://raw.githubusercontent.com/deepchem/deepchem/master/examples/tutorials/assets/atomic_contributions_tutorial_data/logBB_test_.sdf',
    current_dir,
    'logBB_test_.sdf'
)
TEST_DATASET_FILE = os.path.join(current_dir, 'logBB_test_.sdf')
loader = dc.data.SDFLoader(tasks=["p_np"], sanitize=True,
                        featurizer=dc.feat.ConvMolFeaturizer())
test_dataset = loader.create_dataset(TEST_DATASET_FILE, shard_size=2000)
pred =  m.predict(test_dataset)
pred = np.argmax(np.squeeze(pred),axis=1)
ba = metrics.balanced_accuracy_score(y_true=test_dataset.y, y_pred=pred)
print(ba)

平衡精度足够高。现在让我们继续进行模型解释并估计单个原子对预测的贡献。

片段数据集#

现在让我们基于训练集准备一个片段数据集。(也可以使用其他感兴趣的未见过的数据集)。这些片段将被用来评估单个原子的贡献。

对于每个分子，我们将生成一个 ConvMol 对象列表。指定 per_atom_fragmentation=True 告诉它遍历所有的重原子，并特征化当前重原子被删除的分子。也就是说如果某个分子有n个重原子，经过下面过程后，会有n个片段数据，每个片段数据是这个分子中删除了某个原子之后的片段，每个片段数据各不相同。

loader = dc.data.SDFLoader(tasks=[],# dont need task (moreover, passing the task can lead to inconsitencies in data shapes)
                        featurizer=dc.feat.ConvMolFeaturizer(per_atom_fragmentation=True),
                        sanitize=True)
frag_dataset = loader.create_dataset(DATASET_FILE, shard_size=5000)

该数据集仍然具有与原始训练集相同数量的样本，但每个样本现在表示为一个ConvMol对象（每个片段一个）列表，而不是单个ConvMol。

重要提示：片段的顺序取决于输入格式。如果是SDF，则片段顺序与相应 mol 块中的原子顺序相同。如果是SMILES（比如 csv，其中分子用 SMILES 表示），则顺序由 RDKit CanonicalRankAtoms 给出。

print(frag_dataset.X.shape)

我们真的想把每个片段当作一个单独的样本。我们可以使用 FlatteningTransformer 来扁平化片段列表。

tr = dc.trans.FlatteningTransformer(frag_dataset)
frag_dataset = tr.transform(frag_dataset)
print(frag_dataset.X.shape)

预测原子对活性的贡献#

现在我们将预测分子及其片段的活性。然后，对于每个片段，我们将发现活性差异：移除一个原子时活性的变化。

注意：这里，在分类背景下，我们使用模型的概率输出作为活性。所以贡献是概率差，也就是。一个给定的原子增加/减少分子活性的概率。

# whole  molecules
pred = np.squeeze(m.predict(dataset))[:, 1] # probabilitiy of class 1
pred = pd.DataFrame(pred, index=dataset.ids, columns=["Molecule"])  # turn to dataframe for convinience

# fragments
pred_frags = np.squeeze(m.predict(frag_dataset))[:, 1]
pred_frags = pd.DataFrame(pred_frags, index=frag_dataset.ids, columns=["Fragment"])

我们取差值来找出原子的贡献。

# merge 2 dataframes by molecule names
df = pd.merge(pred_frags, pred, right_index=True, left_index=True)
# find contribs
df['Contrib'] = df["Molecule"] - df["Fragment"]

print(df)

我们可以使用 RDKit 的 SimilarityMaps 功能将结果可视化。每个原子都根据它对活性的影响而着色。

def vis_contribs(mols, df, smi_or_sdf = "sdf"):
    # input format of file, which was used to create dataset determines the order of atoms,
    # so we take it into account for correct mapping!
    maps = []
    for mol  in mols:
        wt = {}
        if smi_or_sdf == "smi":
            for n,atom in enumerate(Chem.rdmolfiles.CanonicalRankAtoms(mol)):
                wt[atom] = df.loc[mol.GetProp("_Name"),"Contrib"][n]
        if smi_or_sdf == "sdf":
            for n,atom in enumerate(range(mol.GetNumHeavyAtoms())):
                wt[atom] = df.loc[Chem.MolToSmiles(mol),"Contrib"][n]
        maps.append(SimilarityMaps.GetSimilarityMapFromWeights(mol,wt))
    return maps

让我们看一些图片：

np.random.seed(2000)
maps = vis_contribs(np.random.choice(np.array(mols),10), df)

我们可以看到，芳香族或脂肪族化合物对血脑屏障通透性有积极影响，而极性或带电杂原子对血脑屏障通透性有消极影响。这与文献资料基本一致。

回归的任务#

上面的例子使用了一个分类模型。同样的技术也可以用于回归模型。让我们看一个回归任务，水生毒性（对水生生物梨形虫）。

毒性定义为 log10(IGC50)（抑制菌落生长50%的浓度）。梨形螺旋杆菌的毒性团将通过原子贡献来鉴定。

以上所有步骤都是相同的：加载数据，特征化，构建模型，创建片段数据集，计算贡献，并将它们可视化。

注意：这一次由于是回归，贡献将以活性单位，而不是概率。

current_dir = os.path.dirname(os.path.realpath('__file__'))
dc.utils.download_url(
    'https://raw.githubusercontent.com/deepchem/deepchem/master/examples/tutorials/assets/atomic_contributions_tutorial_data/Tetrahymena_pyriformis_Work_set_OCHEM.sdf',
    current_dir,
    'Tetrahymena_pyriformis_Work_set_OCHEM.sdf'
)
DATASET_FILE =os.path.join(current_dir, 'Tetrahymena_pyriformis_Work_set_OCHEM.sdf')

# create RDKit mol objects, we will need them later
mols = [m for m in Chem.SDMolSupplier(DATASET_FILE) if m is not None ]
loader = dc.data.SDFLoader(tasks=["IGC50"],
                        featurizer=dc.feat.ConvMolFeaturizer(), sanitize=True)
dataset = loader.create_dataset(DATASET_FILE, shard_size=5000)

创建并训练模型。

np.random.seed(2020)
tf.random.set_seed(2020)
m = dc.models.GraphConvModel(1, mode="regression", batch_normalize=False)
m.fit(dataset, nb_epoch=40)

加载测试数据集并检验模型的性能。

current_dir = os.path.dirname(os.path.realpath('__file__'))
dc.utils.download_url(
    'https://raw.githubusercontent.com/deepchem/deepchem/master/examples/tutorials/assets/atomic_contributions_tutorial_data/Tetrahymena_pyriformis_Test_set_OCHEM.sdf',
    current_dir,
    'Tetrahymena_pyriformis_Test_set_OCHEM.sdf'
)




TEST_DATASET_FILE = os.path.join(current_dir, 'Tetrahymena_pyriformis_Test_set_OCHEM.sdf')
loader = dc.data.SDFLoader(tasks=["IGC50"], sanitize= True,
                        featurizer=dc.feat.ConvMolFeaturizer())
test_dataset = loader.create_dataset(TEST_DATASET_FILE, shard_size=2000)
pred = m.predict(test_dataset)
mse = metrics.mean_squared_error(y_true=test_dataset.y, y_pred=pred)
r2 = metrics.r2_score(y_true=test_dataset.y, y_pred=pred)
print(mse)
print(r2)

再次加载训练集，但这次设置为 per_atom_fragmentation=True。

loader = dc.data.SDFLoader(tasks=[], # dont need any task
                        sanitize=True,
                        featurizer=dc.feat.ConvMolFeaturizer(per_atom_fragmentation=True))
frag_dataset = loader.create_dataset(DATASET_FILE, shard_size=5000)
tr = dc.trans.FlatteningTransformer(frag_dataset) # flatten dataset and add ids to each fragment
frag_dataset = tr.transform(frag_dataset)

计算活性差异。

# whole molecules
pred = m.predict(dataset)
pred = pd.DataFrame(pred, index=dataset.ids, columns=["Molecule"])  # turn to dataframe for convenience

# fragments
pred_frags = m.predict(frag_dataset)
pred_frags = pd.DataFrame(pred_frags, index=frag_dataset.ids, columns=["Fragment"])  # turn to dataframe for convenience

# merge 2 dataframes by molecule names
df = pd.merge(pred_frags, pred, right_index=True, left_index=True)
# find contribs
df['Contrib'] = df["Molecule"] - df["Fragment"]

让我们取一些活性中等（不是非常活跃/不活跃）的分子，并可视化其原子贡献。

maps = vis_contribs([mol for mol in mols if float(mol.GetProp("IGC50"))>3 and float(mol.GetProp("IGC50"))<4][:10], df)

我们可以看到，已知的毒株为绿色，分别为硝基芳烃类、卤代芳烃类、长烷基链类和醛类：而羧基、醇和氨基则具有解毒作用，与文献[3]一致。

附录#

在本教程中，我们操作的是 SDF 文件。然而，如果我们使用 SMILES 作为输入的 CSV 文件，则 dataframe 中原子的顺序与原始原子的顺序不对应。如果我们想要恢复每个分子的原始原子顺序（为了添加到我们的 dataframe 中），我们需要使用 RDKit 的 Chem.rdmolfiles.CanonicalRankAtoms。这里有一些实用程序可以做到这一点。

我们可以添加一个包含原子 id 的列（就像输入分子那样），并将得到的 dataframe 用于“python-rdkit-deepchem”环境之外的任何其他软件的分析。

def get_mapping(mols, mol_names):
    """perform mapping:
    atom number original <-> atom number(position)
    after ranking (both 1-based)"""
    # mols - RDKit mols
    # names  - any seq of strings
    # return list of nested lists: [[molecule, [atom , atom, ..], [...]]
    assert(len(mols)==len(mol_names))
    mapping = []
    for m,n in zip(mols, mol_names):
        atom_ids = [i+1 for i in list(Chem.rdmolfiles.CanonicalRankAtoms(m))]
        mapping.append([n, atom_ids])
    return mapping

def append_atomid_col(df, mapping):
    # add column with CORRECT atom number(position)
    for i in mapping:
        df.loc[i[0],"AtomID"] = i[1]
    return df

参考文献：#

Polishchuk, P., O. Tinkov, T. Khristova, L. Ognichenko, A. Kosinskaya, A. Varnek & V. Kuz’min (2016) Structural and Physico-Chemical Interpretation (SPCI) of QSAR Models and Its Comparison with Matched Molecular Pair Analysis. Journal of Chemical Information and Modeling, 56, 1455-1469.
Riniker, S. & G. Landrum (2013) Similarity maps - a visualization strategy for molecular fingerprints and machine-learning methods. Journal of Cheminformatics, 5, 43.
1. Matveieva, M. T. D. Cronin, P. Polishchuk, Mol. Inf. 2019, 38, 1800084.
Matveieva, M., Polishchuk, P. Benchmarks for interpretation of QSAR models. J Cheminform 13, 41 (2021). https://doi.org/10.1186/s13321-021-00519-x

使用 Trident Chemwidgets 进行交互式模型评估#