分子表征、自监督与蛋白-配体建模

编者导语

本板块覆盖分子机器学习的核心技术路线：分子如何被编码成机器可用的表示 （SMILES、图结构、指纹与描述符）、神经网络如何处理这些表示、 自监督与预训练如何缓解标注稀缺、以及蛋白-配体相互作用的建模。 板块含 149 张图，是全书图量最集中的部分之一；公式不多（7 个）， 因为方法叙述远多于数学推导。 几处需要读者自己设限的地方：自监督模型的收益高度依赖下游任务的选择、 迁移并非自动发生；以同一复合物的不同姿态划分训练与测试会显著高估模型能力； 计算亲和力预测不能替代实验测定。 本书对这一整编的定位是：深度学习的成功案例集中在特定任务上， 把它当通用替代方案的期望并不现实。
前置：ch51。

Creating Models with TensorFlow and PyTorch

In the tutorials so far, we have used standard models provided by DeepChem. This is fine for many applications, but sooner or later you will want to create an entirely new model with an architecture you define yourself. DeepChem provides integration with both TensorFlow (Keras) and PyTorch, so you can use it with models from either of these frameworks.

Colab

This tutorial and the rest in this sequence are designed to be done in Google colab. If you’d like to open this notebook in colab, you can use the following link.

1

In [ ]: !pip install --pre deepchem

There are actually two different approaches you can take to using TensorFlow or PyTorch models with DeepChem. It depends on whether you want to use TensorFlow/PyTorch APIs or DeepChem APIs for training and evaluating your model. For the former case, DeepChem’s Dataset class has methods for easily adapting it to use with other frameworks. make_tf_dataset() returns a tensorflow.data.Dataset object that iterates over the data. make_pytorch_dataset() returns a torch.utils.data.IterableDataset that iterates over the data. This lets you use DeepChem’s datasets, loaders, featurizers, transformers, splitters, etc. and easily integrate them into your existing TensorFlow or PyTorch code.

But DeepChem also provides many other useful features. The other approach, which lets you use those features, is to wrap your model in a DeepChem Model object. Let’s look at how to do that.

KerasModel

KerasModel is a subclass of DeepChem’s Model class. It acts as a wrapper around a tensorflow.keras.Model . Let’s see an example of using it. For this example, we create a simple sequential model consisting of two dense layers.

1
2
3
4
5
6
7
8
9
10

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

For this example, we used the Keras Sequential class. Our model consists of a dense layer with ReLU activation, 50% dropout to provide regularization, and a final layer that produces a scalar output. We also need to specify the loss function to use when training the model, in this case L₂ loss. We can now train and evaluate the model exactly as we would with any other DeepChem model. For example, let’s load the Delaney solubility dataset. How does our model do at predicting the solubilities of molecules based on their extended-connectivity fingerprints (ECFPs)?

1
2
3
4
5
6
7
8
9
10

In [ ]:
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]))

training set score: {'pearson_r2_score': 0.9766804253639305}
test set score: {'pearson_r2_score': 0.7048814451615332}

TorchModel

TorchModel works just like KerasModel , except it wraps a torch.nn.Module . Let’s use PyTorch to create another model just like the previous one and train it on the same data.

1
2
3
4

import torch
pytorch_model = torch.nn.Sequential(
    torch.nn.Linear(1024, 1000),
    torch.nn.ReLU(),

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

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)

1
2
3
4
5
6
7
8
9
10

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

training set score: {'pearson_r2_score': 0.9760781898204121}
test set score: {'pearson_r2_score': 0.6981331812360332}

Computing Losses

Now let’s see a more advanced example. In the above models, the loss was computed directly from the model’s output. Often that is fine, but not always. Consider a classification model that outputs a probability distribution. While it is possible to compute the loss from the probabilities, it is more numerically stable to compute it from the logits.

To do this, we create a model that returns multiple outputs, both probabilities and logits. KerasModel and TorchModel let you specify a list of “output types”. If a particular output has type ‘prediction’ , that means it is a normal output that should be returned when you call predict() . If it has type ‘loss’ , that means it should be passed to the loss function in place of the normal outputs.

Sequential models do not allow multiple outputs, so instead we use a subclassing style model.

We can train our model on the BACE dataset. This is a binary classification task that tries to predict whether a molecule will inhibit the enzyme BACE-1.

1
2
3
4
5
6

In [ ]: 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]))

1
2

training set score: {'roc_auc_score': 0.9995809177900399}
test set score: {'roc_auc_score': 0.7629528985507246}

Similarly, we will create a custom Classifier Model class to be used with TorchModel . Using similar reasoning to the above KerasModel , a custom model allows for easy capturing of the unscaled output (logits in Tensorflow) of the second dense layer. The custom class allows definition of how forward pass is done; enabling capture of the logits right before the final sigmoid is applied to produce the prediction

Finally, an instance of ClassificationModel is coupled with a loss function that requires both the prediction and logits to produce an instance of TorchModel to train.

1
2
3
4
5
6
7
8
9
10
11
12

In [ ] : 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

1
2
3

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

We will use the same BACE dataset. As before, the model will try to do a binary classification task that tries to predict whether a molecule will inhibit the enzyme BACE-1.

1
2
3
4
5
6

In [ ]: 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]))

1
2

training set score: {'roc_auc_score': 0.9996340015366347}
test set score: {'roc_auc_score': 0.7615036231884058}

Other Features

KerasModel and TorchModel have lots of other features. Here are some of the more important ones.

Automatically saving checkpoints during training.

Logging progress to the console, to TensorBoard, or to Weights & Biases.

Custom loss functions that you define with a function of the form f(outputs, labels, weights)

Early stopping using the ValidationCallback class.

Loading parameters from pre-trained models.

Estimating uncertainty in model outputs.

Identifying important features through saliency mapping.

By wrapping your own models in a KerasModel or TorchModel , you get immediate access to all these features. See the API documentation for full details on them.

Congratulations! Time to join the Community!

Congratulations on completing this tutorial notebook! If you enjoyed working through the tutorial, and want to continue working with DeepChem, we encourage you to finish the rest of the tutorials in this series. You can also help the DeepChem community in the following ways:

Star DeepChem on GitHub

This helps build awareness of the DeepChem project and the tools for open source drug discovery that we’re trying to build.

Join the DeepChem Discord

The DeepChem Discord hosts a number of scientists, developers, and enthusiasts interested in deep learning for the life sciences. Join the conversation!

Citing This Tutorial

If you found this tutorial useful please consider citing it using the provided BibTeX.

1
2
3
4
5
6

@manual{Intro1,
    title={5},
    organization={DeepChem},
    author={Ramsundar, Bharath and Rebel, Alles},
    howpublished = {\url{https://github.com/deepchem/deepchem/blob/master/examples/tutorials/Creating_Models_with_year={2021},
}

Introduction to Graph Convolutions

In this tutorial we will learn more about “graph convolutions.” These are one of the most powerful deep learning tools for working with molecular data. The reason for this is that molecules can be naturally viewed as graphs.

Note how standard chemical diagrams of the sort we’re used to from high school lend themselves naturally to visualizing molecules as graphs. In the remainder of this tutorial, we’ll dig into this relationship in significantly more detail. This will let us get a deeper understanding of how these systems work.

Colab

This tutorial and the rest in this sequence are designed to be done in Google colab. If you’d like to open this notebook in colab, you can use the following link.

Open in Colab

!pip install –pre deepchem

What are Graph Convolutions?

Consider a standard convolutional neural network (CNN) of the sort commonly used to process images. The input is a grid of pixels. There is a vector of data values for each pixel, for example the red, green, and blue color channels. The data passes through a series of convolutional layers. Each layer combines the data from a pixel and its neighbors to produce a new data vector for the pixel. Early layers detect small scale local patterns, while later layers detect larger, more abstract patterns. Often the convolutional layers alternate with pooling layers that perform some operation such as max or min over local regions.

Graph convolutions are similar, but they operate on a graph. They begin with a data vector for each node of the graph (for example, the chemical properties of the atom that node represents). Convolutional and pooling layers combine information from connected nodes (for example, atoms that are bonded to each other) to produce a new data vector for each node.

Training a GraphConvModel

Let’s use the MoleculeNet suite to load the Tox21 dataset. To featurize the data in a way that graph convolutional networks can use, we set the featurizer option to ‘GraphConv’ . The MoleculeNet call returns a training set, a validation set, and a test set for us to use. It also returns tasks , a list of the task names, and transformers , a list of data transformations that were applied to preprocess the dataset. (Most deep networks are quite finicky and require a set of data transformations to ensure that training proceeds stably.)

Note: While importing deepchem, if you see any warnings, ignore them for now. Deepchem is a vast library and there are many things that can cause minor warnings to occur. Almost always, it doesn’t require any action from your side.

1
2
3

In [4]: import deepchem as dc
tasks, datasets, transformers = dc.molnet.load_tox21(featurizer='GraphConv')
train_dataset, valid_dataset, test_dataset = datasets

Let’s now train a graph convolutional network on this dataset. DeepChem has the class GraphConvModel that wraps a standard graph convolutional architecture underneath the hood for user convenience. Let’s instantiate an object of this class and train it on our dataset.

1
2
3

In [5]: n_tasks = len(tasks)
num_features = train_dataset.X[0].get_atom_features().shape[1]
model = dc.models.torch_models.GraphConvModel(n_tasks, mode='classification', number_input_features=[num_features], model.fit(train_dataset, nb_epoch=50)

1

Out[5]: 0.29102970123291017

Let’s try to evaluate the performance of the model we’ve trained. For this, we need to define a metric, a measure of model performance. dc.metrics holds a collection of metrics already. For this dataset, it is standard to use the ROC AUC score, the area under the receiver operating characteristic curve (which measures the tradeoff between precision and recall). Luckily, the ROC-AUC score is already available in DeepChem.

To measure the performance of the model under this metric, we can use the convenience function model.evaluate() .

1
2
3
4
5
6

In [6]: 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))

Training set score: {'roc_auc_score': 0.970785822904073}
Test set score: {'roc_auc_score': 0.7112009940440461}

The results are pretty good, and GraphConvModel is very easy to use. But what’s going on under the hood? Could we build GraphConvModel ourselves? Of course! DeepChem provides Keras layers for all the calculations involved in a graph convolution. We are going to apply the following layers from DeepChem.

GraphConv layer: This layer implements the graph convolution. The graph convolution combines per-node feature vectures in a nonlinear fashion with the feature vectors for neighboring nodes. This “blends” information in local neighborhoods of a graph.

GraphPool layer: This layer does a max-pooling over the feature vectors of atoms in a neighborhood. You can think of this layer as analogous to a max-pooling layer for 2D convolutions but which operates on graphs instead.

GraphGather : Many graph convolutional networks manipulate feature vectors per graph-node. For a molecule for example, each node might represent an atom, and the network would manipulate atomic feature vectors that summarize the local chemistry of the atom. However, at the end of the application, we will likely want to work with a molecule level feature representation. This layer creates a graph level feature vector by combining all the node level feature vectors.

Apart from this we are going to apply standard neural network layers such as Dense, BatchNormalization and Softmax layer.

Training a custom Graph Convolution network

As you may have seen in the previous tutorials, DeepChem offers both PyTorch and Tensorflow functionalities. However, most of our work moving forward will leverage the PyTorch ecosystem.

Let’s look at the Tensorflow implementation first.

Tensorflow

1
2
3
4
5
6
7
8
9
10
11
12

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

batch_size = 100

class GraphConvModelTensorflow(tf.keras.Model):
    def __init__(self):
    super(GraphConvModelTensorflow, self).__init__()
    self.gc1 = GraphConv(128, activation_fn=tf.nn.tanh)
    self.batch_norm1 = layers.BatchNorm()
    self.gp1 = GraphPool()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

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

self.dense1 = layers.Dense(256, activation=tf.nn.tanh)
self.batch_norm3 = layers.BatchNorm()
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-unit])
    dense1_output = self.dense1(gp2_output)
    batch_norm3_output = self.batch_norm3(dense1_output)
    readout_output = self.readout([batch_norm3_output] + inputs[1-unit])
    logits_output = self.logits(self.dense2(readout_output))
    return self.softmax(logits_output)

We can now see more clearly what is happening. There are two convolutional blocks, each consisting of a GraphConv followed by batch normalization, followed by a GraphPool to do max pooling. We finish up with a dense layer, another batch normalization, a GraphGather to combine the data from all the different nodes, and a final dense layer to produce the global output.

Let’s now create the DeepChem model which will be a wrapper around the Keras model that we just created. We will also specify the loss function so the model knows the objective to minimize.

1

In [8]: model = dc.models.KerasModel(GraphConvModelTensorflow(), loss=dc.models.losses.CategoricalCrossEntropy())

What are the inputs to this model? A graph convolution requires a complete description of each molecule, including the list of nodes (atoms) and a description of which ones are bonded to each other. In fact, if we inspect the dataset we see that the feature array contains Python objects of type ConvMol .

1

In [9]: test_dataset.X[0]

1

Out[9]: <deepchem.feat.mol_graphs.ConvMol at 0x7bf66bfa1160>

Models expect arrays of numbers as their inputs, not Python objects. We must convert the ConvMol objects into the particular set of arrays expected by the GraphConv , GraphPool , and GraphGather layers. Fortunately, the ConvMol class includes the code to do this, as well as to combine all the molecules in a batch to create a single set of arrays.

The following code creates a Python generator that given a batch of data generates the lists of inputs, labels, and weights whose values are Numpy arrays. atom_features holds a feature vector of length 75 for each atom. The other inputs are required to support minibatching in TensorFlow. degree_slice is an indexing convenience that makes it easy to locate atoms from all molecules with a given degree. membership determines the membership of atoms in molecules (atom i belongs to molecule membership[i] ). deg_adjs is a list that contains adjacency lists grouped by atom degree. For more details, check out the code.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

In [10]:
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)

Now, we can train the model using fit_generator(generator) which will use the generator we’ve defined to train the model.

1

In [11]: model.fit_generator(data_generator(train_dataset, epochs=50))

1

Number of input features: 75

1

Out[11]: 0.23354644775390626

Now that we have trained our graph convolutional method, let’s evaluate its performance. We again have to use our defined generator to evaluate model performance.

1
2

In [12]: 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))

PyTorch

Before working on the PyTorch implementation, we must import a few crucial layers from the torch_models collection. These are PyTorch implementations of GraphConv , GraphPool and GraphGather which we used in the tensorflow’s implementation as well.

1
2
3

In [13]: import torch
import torch.nn as nn
from deepchem.models.torch_models.layers import GraphConv, GraphGather, GraphPool

PyTorch’s GraphConv requires the number of input features to be specified, hence we can extract that piece of information by following steps:

1. First we get a sample from the dataset.

2. Next we slice and separate the node_features (which is the first element of the list, hence the index 0).

3. Finally, we obtain the number of features by finding the shape of the array.

1
2
3
4

In [14]: sample_batch = next(data_generator(train_dataset))
node_features = sample_batch[0][0]
num_input_features = node_features.shape[1]
print(f"Number of input features: {num_input_features}")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

class GraphConvModelTorch(nn.Module):
    def __init__(self):
    super(GraphConvModelTorch, self).__init__()

    self.gc1 = GraphConv(out_channel=128, number_input_features=num_input_features, activation_fn=nn.Tanh()
    self.batch_norm1 = nn.BatchNorm1d(128)
    self.gp1 = GraphPool()

    self.gc2 = GraphConv(out_channel=128, number_input_features=128, activation_fn=nn.Tanh())
    self.batch_norm2 = nn.BatchNorm1d(128)
    self.gp2 = GraphPool()

    self.dense1 = nn.Linear(128, 256)
    self.act3 = nn.Tanh()
    self.batch_norm3 = nn.BatchNorm1d(256)
    self.readout = GraphGather(batch_size=batch_size, activation_fn=nn.Tanh())

    self.dense2 = nn.Linear(512, n_tasks * 2)

    self.logits = lambda data: data.view(-1, n_tasks, 2)
    self.softmax = nn.softmax(dim=-1)

def forward(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_norm2(gc2_output)
    gp2_output = self.gp2([batch_norm2_output] + inputs[1-unit])
    dense1_output = self.act3(self.dense1(gp2_output))
    batch_norm3_output = self.batch_norm3(dense1_output)
    readout_output = self.readout([batch_norm3_output] + inputs[1-unit])

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

1

In [18]: model = dc.models.TorchModel(GraphConvModelTorch(), loss=dc.models.losses.CategoricalCrossEntropy()) 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))

1
2

Training set score: {'roc_auc_score': 0.9838238607233897}
Test set score: {'roc_auc_score': 0.6923516284964811}

Success! Both the models we’ve constructed behave nearly identically to GraphConvModel . If you’re looking to build your own custom models, you can follow the examples we’ve provided here to do so. We hope to see exciting constructions from your end soon!

Congratulations! Time to join the Community!

Congratulations on completing this tutorial notebook! If you enjoyed working through the tutorial, and want to continue working with DeepChem, we encourage you to finish the rest of the tutorials in this series. You can also help the DeepChem community in the following ways:

Star DeepChem on GitHub

This helps build awareness of the DeepChem project and the tools for open source drug discovery that we’re trying to build.

Join the DeepChem Gitter

The DeepChem Gitter hosts a number of scientists, developers, and enthusiasts interested in deep learning for the life sciences. Join the conversation!

1
2
3
4
5
6
7

In [1]: !pip install -qq --pre deepchem
import deepchem
import warnings
warnings.filterwarnings('ignore')
deepchem.__version__

Out[1]: '2.6.0.dev'

Going Deeper On Molecular Featurizations

One of the most important steps of doing machine learning on molecular data is transforming the data into a form amenable to the application of learning algorithms. This process is broadly called “featurization” and involves turning a molecule into a vector or tensor of some sort. There are a number of different ways of doing that, and the choice of featurization is often dependent on the problem at hand. We have already seen two such methods: molecular fingerprints, and ConvMol objects for use with graph convolutions. In this tutorial we will look at some of the others.

Colab

This tutorial and the rest in this sequence can be done in Google colab. If you’d like to open this notebook in colab, you can use the following link.

Featurizers

In DeepChem, a method of featurizing a molecule (or any other sort of input) is defined by a Featurizer object. There are three different ways of using featurizers.

1. When using the MoleculeNet loader functions, you simply pass the name of the featurization method to use. We have seen examples of this in earlier tutorials, such as featurizer=’ECFP’ or featurizer=’GraphConv’

2. You also can create a Featurizer and directly apply it to molecules. For example:

1
2
3
4
5
6
7
8
9
10

In [2]: import deepchem as dc
    featurizer = dc.feat.CircleFingerprint()
    print(featurizer(['CC', 'CCC', 'CCO']))

[[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]]
[09:05:59] DEPRECATION WARNING: please use MorganGenerator
[09:05:59] DEPRECATION WARNING: please use MorganGenerator
[09:05:59] DEPRECATION WARNING: please use MorganGenerator

3. When creating a new dataset with the DataLoader framework, you can specify a Featurizer to use for processing the data. We will see this in a future tutorial.

We use propane (CH₃CH₂CH₃, represented by the SMILES string ‘CCC’ ) as a running example throughout this tutorial. Many of the featurization methods use conformers of the molecules. A conformer can be generated using the ConformerGenerator class in deepchem.utils.conformers .

RDKitDescriptors

RDKitDescriptors featurizes a molecule by using RDKit to compute values for a list of descriptors. These are basic physical and chemical properties: molecular weight, polar surface area, numbers of hydrogen bond donors and acceptors, etc. This is most useful for predicting things that depend on these high level properties rather than on detailed molecular structure.

Intrinsic to the featurizer is a set of allowed descriptors, which can be accessed using RDKitDescriptors.allowedDescriptors . The featurizer uses the descriptors in rdkit.Chem.Descriptors.descList , checks if they are in the list of allowed descriptors, and computes the descriptor value for the molecule.

Let’s print the values of the first ten descriptors for propane.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

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

MaxAbsEStateIndex 2.125
MaxEStateIndex 2.125
MinAbsEStateIndex 1.25
MinEStateIndex 1.25
qed 0.3854706587740357
SPS 6.0
MolWt 44.097
HeavyAtomMolWt 36.033
ExactMolWt 44.062600255999996
NumValenceElectrons 20.0

[09:07:17] DEPRECATION WARNING: please use MorganGenerator
[09:07:17] DEPRECATION WARNING: please use MorganGenerator
[09:07:17] DEPRECATION WARNING: please use MorganGenerator

Of course, there are many more descriptors than this.

1

In [4]: print('The number of descriptors present is: ', len(features))

1

The number of descriptors present is: 210

WeaveFeaturizer and MolGraphConvFeaturizer

We previously looked at graph convolutions, which use ConvMolFeaturizer to convert molecules into ConvMol objects. Graph convolutions are a special case of a large class of architectures that represent molecules as graphs. They work in similar ways but vary in the details. For example, they may associate data vectors with the atoms, the bonds connecting them, or both. They may use a variety of techniques to calculate new data vectors from those in the previous layer, and a variety of techniques to compute molecule level properties at the end.

DeepChem supports lots of different graph based models. Some of them require molecules to be featurized in slightly different ways. Because of this, there are two other featurizers called WeaveFeaturizer and MolGraphConvFeaturizer . They each convert molecules into a different type of Python object that is used by particular models. When using any graph based model, just check the documentation to see what featurizer you need to use with it.

CoulombMatrix

All the models we have looked at so far consider only the intrinsic properties of a molecule: the list of atoms that compose it and the bonds connecting them. When working with flexible molecules, you may also want to consider the different conformations the molecule can take on. For example, when a drug molecule binds to a protein, the strength of the binding depends on specific interactions between pairs of atoms. To predict binding strength, you probably want to consider a variety of possible conformations and use a model that takes them into account when making predictions.

The Coulomb matrix is one popular featurization for molecular conformations. Recall that the electrostatic Coulomb interaction between two charges is proportional to

matrix where each element gives the strength of the electrostatic interaction between two atoms. It contains information both about the charges on the atoms and the distances between them. More information on the functiona forms used can be found here.

To apply this featurizer, we first need a set of conformations for the molecule. We can use the ConformerGenerator class to do this. It takes a RDKit molecule, generates a set of energy minimized conformers, and prunes the set to only include ones that are significantly different from each other. Let’s try running it for propane.

1
2
3
4
5

In [5]: 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()))

1
2
3
4

In [6]: butane_mol = generator.generate_conformers(Chem.MolFromSmiles('CCCC'))
print("Number of available conformers for butane: ", len(butane_mol.GetConformers()))
Number of available conformers for butane: 3
Now we can create a Coulomb matrix for our molecule.

Number of available conformers for propane: 1

It only found a single conformer. This shouldn’t be surprising, since propane is a very small molecule with hardly any flexibility. Let’s try adding another carbon.

1
2
3

In [7]: coulomb_mat = dc.feat.CoulombMatrix(max_atoms=20)
features = coulomb_mat(propane_mol)
print(features)

Notice that many elements are 0. To combine multiple molecules in a batch we need all the Coulomb matrices to be the same size, even if the molecules have different numbers of atoms. We specified max_atoms=20 , so the returned matrix

has size (20, 20). The molecule only has 11 atoms, so only an 11 by 11 submatrix is nonzero.

CoulombMatrixEig

An important feature of Coulomb matrices is that they are invariant to molecular rotation and translation, since the interatomic distances and atomic numbers do not change. Respecting symmetries like this makes learning easier. Rotating a molecule does not change its physical properties. If the featurization does change, then the model is forced to learn that rotations are not important, but if the featurization is invariant then the model gets this property automatically.

Coulomb matrices are not invariant under another important symmetry: permutations of the atoms’ indices. A molecule’s physical properties do not depend on which atom we call “atom 1”, but the Coulomb matrix does. To deal with this, the CoulumbMatrixEig featurizer was introduced, which uses the eigenvalue spectrum of the Coulumb matrix and is invariant to random permutations of the atom’s indices. The disadvantage of this featurization is that it contains much less information (

N
eigenvalues instead of an
N \times N
matrix), so models will be more limited in what they can learn.

CoulombMatrixEig inherits from CoulombMatrix and featurizes a molecule by first computing the Coulomb matrices for different conformers of the molecule and then computing the eigenvalues for each Coulomb matrix. These eigenvalues are then padded to account for variation in number of atoms across molecules.

1
2
3
4
5
6
7
8

In [8]: coulomb_mat_eig = dc.feat.CoulombMatrixEig(max_atoms=20)
features = coulomb_mat_eig(propane_mol)
print(features)

[[60.07620303 29.62963149 22.75497781 0.5713786 0.28781332 0.28548338
0.27558187 0.18163794 0.17460999 0.17059719 0.16640098 0.
0.    0.    0.    0.    0.    0.
0.    0.    ]]

SMILES Tokenization and Numericalization

So far, we have looked at featurization techniques that translate the implicit structural and physical information in SMILES data into more explicit features that help our models learn and make predictions. In this section, we will preprocess SMILES strings into a format that can be fed to sequence models, such as 1D convolutional neural networks and transformers, and enables these models to learn their own representations of molecular properties.

To prepare SMILES strings for a sequence model, we break them down into lists of substrings (called tokens) and turn them into lists of integer values (numericalization). Sequence models use those integer values as indices of an embedding matrix, which contains a vector of floating-point numbers for each token in the vocabulary. These embedding vectors are updated during model training. This process allows the sequence model to learn its own representations of the molecular properties implicit in the training data.

We will use DeepChem’s BasicSmilesTokenizer and the Tox21 dataset from MoleculeNet to demonstrate the process of tokenizing SMILES.

1
2
3
4
5
6
7
8
9

In [9]: import numpy as np

In [10]: tasks, datasets, transformers = dc.molnet.load_tox21(featurizer="Raw")
train_dataset, valid_dataset, test_dataset = datasets
print(train_dataset)

<DiskDataset X.shape: (6264,), y.shape: (6264, 12), w.shape: (6264, 12), task_names: ['NR-AR' 'NR-AR-LBD' 'NR-Ah R' ... 'SR-HSE' 'SR-MMP' 'SR-p53']>

We loaded the datasets with featurizer="Raw". Now we obtain the SMILES from their ids attributes.

1
2
3
4
5
6
7
8
9
10

In [11]: train_smiles = train_dataset.ids
valid_smiles = valid_dataset.ids
test_smiles = test_dataset.ids
print(train_smiles[:5])

['CC(0)(P(=0)(0)0)P(=0)(0)0' 'CC(C)(C)00C(C)(C)CCC(C)(C)00C(C)(C)C'
'0C[C@H](0)[C@@H](0)[C@H](0)C0'
'CCCCCC(=0)[0-].CCCCCC(=0)[0-].[Zn+2]' 'CC(C)COC(=0)C(C)C']

Next we define our tokenizer and map it onto all our data to convert the SMILES strings into lists of tokens. The BasicSmileTokenizer breaks down SMILES roughly at atom level.

1

In [12]: tokenizer = dc.feat.smiles_tokenizer.BasicSmilesTokenizer()

1

[19, 19, 2, 24, 3, 2, 25, 2, 16, 24, 3, 2, 24, 3, 24, 3, 25, 2, 16, 24, 3, 2, 24, 3, 24]

1
2
3
4
5
6
7

train_tok = list(map(tokenizer.tokenize, train_smiles))
valid_tok = list(map(tokenizer.tokenize, valid_smiles))
test_tok = list(map(tokenizer.tokenize, test_smiles))
print(train_tok[0])
len(train_tok)

['C', 'C', '(', '0', ')' ', '(', 'P', '(', '=', '0', ')' ', '(', '0', ')' ', '0', ')' ', 'P', '(', '=', '0', ')' ', '(', '0', ')' ', '0']

1

Out[12]: 6264

Now we have tokenized versions of all SMILES strings in our dataset. To convert those into lists of integer values we first need to create a list of all possible tokens in our dataset. That list is called the vocabulary. We also add the empty string “” to our vocabulary in order to correctly handle trailing zeros when decoding zero-padded numericalized SMILES.

1
2
3
4
5
6
7

In [13]: flatten = lambda l: [item for items in l for item in items]
    all_toks = flatten(train_tok) + flatten(valid_tok) + flatten(test_tok)
    vocab = sorted(set(all_toks + [""]))
    print(vocab[:12], "...", vocab[-12:])
    len(vocab)

['', '#', '(', ')', '-', '.', '/', '1', '2', '3', '4', '5'] ... ['[n+]', '[n-]', '[nH+]', '[nH]', '[o+]', '[s+]', '[se]', '\', 'c', 'n', 'o', 's']

1

Out[13]: 128

To numericalize tokenized SMILES strings we create a str2int dictionary which assigns a number to each token in the dictionary. We also create the reverse int2str dictionary and define the corresponding encode and decode functions. Finally we map the encode function on the tokenized data to obtain numericalized SMILES data.

1
2
3
4

In [14]: str2int = {s:i for i, s in enumerate(vocab)}
int2str = {i:s for i, s in enumerate(vocab)}
print(f"str2int: {dict(list(str2int.items())[5])} ...")
print(f"int2str: {dict(list(int2str.items())[5])} ...")

1
2
3
4

In [16]: train_num = list(map(encoder, train_tok))
valid_num = list(maprika, valid_tok))
test_num = list(map(encode, test_tok))
print(train_num[0])

Lastly, we would like to combine all molecules in a dataset in an np.array so they can be served to a model in batches. To achieve that, all sequences have to be of the same length. As in the CoulombMatrix section, we achieve that by appending zeros up to a fixed value.

1
2

In [17]: max_len = max(map(len, train_num + valid_num + test_num))
max_len

1

Out[17]: 240

The longest sequence across all Tox21 datasets has length 240 , so we use that as our fixed length. We create a zero_pad function, map it to all numericalized SMILES, and turn them into np.array s.

1
2
3
4
5
6
7
8
9
10
11
12

In [18]: zero_pad = lambda x: x + [0] * (max_len - len(x))
train_numpad = np.array(list(map(zero_pad, train_num)))
valid_numpad = np.array(list(map(zero_pad, valid_num)))
test_numpad = np.array(list(map(zero_pad, test_num)))
train_numpad
Out[18]: array([[19, 19, 2, ..., 0, 0, 0],
    [19, 19, 2, ..., 0, 0, 0],
    [24, 19, 42, ..., 0, 0, 0],
    ...,
    [24, 16, 19, ..., 0, 0, 0],
    [19, 19, 2, ..., 0, 0, 0],
    [19, 19, 2, ..., 0, 0, 0]])

We can check that the zero-padded data still converts back to the correct SMILES string using the decode function on a random datapoint.

1
2
3
4
5

In [19]: idx = np.random.randint(0, train_numpad.shape[0], 1).item()
print(train_smiles[idx])
print("".join(decode(train_numpad[idx])))
Cclcc(C(C)(C)c2ccc(0)c(C)c2)ccc10
Cclcc(C(C)(C)c2ccc(0)c(C)c2)ccc10

The padded data passes the test. It is now in the correct format to be used for training of a sequence model, but it doesn’t yet interface nicely with DeepChem’s training framework. To change that, we define a tokenize_smiles function that combines all the steps spelled out above to process a single datapoint. Additionally, we define a SmilesFeaturizer that uses our custom tokenize_smiles function in its _featurize method and instanciate it as smiles_featurizer passing it our vocab and max_len .

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

In [20]: def tokenize_smiles(x, vocab, max_len):
    tokenizer = dc.feat.smiles_tokenizer.BasicSmilesTokenizer()
    str2int = {s:i for i, s in enumerate(vocab)}
    encode = lambda s: [str2int[tok] for tok in s]
    zero_pad = lambda x: x + [0] * (max_len - len(x))
    x = tokenizer.tokenize(x)
    x = encode(x)
    x = zero_pad(x)
    return np.array(x)

class SmilesFeaturizer(dc.feat.Featurizer):
    def __init__(self, feat_func, vocab, max_len):
    self.feat_func = feat_func
    self.vocab = vocab
    self.max_len = max_len

    def _featurize(self, x):
    return self.feat_func(x, self.vocab, self.max_len)

smiles_featurizer = SmilesFeaturizer(tokenize_smiles, vocab, max_len)

Finally, we use the smiles_featurizer to create new Tox21 datasets that contain tokenized and numericalized SMILES in their X attribute.

1

In [21]: tasks, datasets, transformers = dc.molnet.load_tox21(featurizer=smiles_featurizer) print(datasets[0].X)

1
2
3
4
5
6
7
8

[09:24:48] WARNING: not removing hydrogen atom without neighbors
[[19 19 2 ... 0 0 0]
[19 19 2 ... 0 0 0]
[24 19 42 ... 0 0 0]
...
[24 16 19 ... 0 0 0]
[19 19 2 ... 0 0 0]
[19 19 2 ... 0 0 0]]

The datasets are now ready to be used with your custom DeepChem sequence model. Don’t forget to wrap your model into the appropriate DeepChem model class.

Congratulations! Time to join the Community!

Congratulations on completing this tutorial notebook! If you enjoyed working through the tutorial, and want to continue working with DeepChem, we encourage you to finish the rest of the tutorials in this series. You can also help the DeepChem community in the following ways:

Star DeepChem on GitHub

This helps build awareness of the DeepChem project and the tools for open source drug discovery that we’re trying to build.

Join the DeepChem Gitter

The DeepChem Gitter hosts a number of scientists, developers, and enthusiasts interested in deep learning for the life sciences. Join the conversation!

Citing This Tutorial

If you found this tutorial useful please consider citing it using the provided BibTeX.

1
2
3
4
5
6
7

In [ ]: @manual{Intro7,
title={Going Deeper on Molecular Featurizations},
organization={DeepChem},
author={Ramsundar, Bharath},
howpublished = {\url{https://github.com/deepchem/deepchem/blob/master/examples/tutorials/Going_Deeper_on_Molec
year={2021},
}

Working With Splitters

When using machine learning, you typically divide your data into training, validation, and test sets. The MoleculeNet loaders do this automatically. But how should you divide up the data? This question seems simple at first, but it turns out to be quite complicated. There are many ways of splitting up data, and which one you choose can have a big impact on the reliability of your results. This tutorial introduces some of the splitting methods provided by DeepChem.

Colab

This tutorial and the rest in this sequence can be done in Google colab. If you’d like to open this notebook in colab, you can use the following link.

Open in Colab

1
2
3

!pip install --pre deepchem
import deepchem
deepchem.__version__

‘2.7.1’

Splitters

In DeepChem, a method of splitting samples into multiple datasets is defined by a Splitter object. Choosing an appropriate method for your data is very important. Otherwise, your trained model may seem to work much better than it really does.

Consider a typical drug development pipeline. You might begin by screening many thousands of molecules to see if they bind to your target of interest. Once you find one that seems to work, you try to optimize it by testing thousands of minor variations on it, looking for one that binds more strongly. Then perhaps you test it in animals and find it has unacceptable toxicity, so you try more variations to fix the problems.

This has an important consequence for chemical datasets: they often include lots of molecules that are very similar to each other. If you split the data into training and test sets in a naive way, the training set will include many molecules that are very similar to the ones in the test set, even if they are not exactly identical. As a result, the model may do very well on the test set, but then fail badly when you try to use it on other data that is less similar to the training data.

Let’s take a look at a few of the splitters found in DeepChem.

• General Splitters

○ RandomSplitter

○ RandomGroupSplitter

○ RandomStratifiedSplitter

○ SingletaskStratifiedSplitter

○ IndexSplitter

○ SpecifiedSplitter

○ TaskSplitter

• Molecular Splitters

○ ScaffoldSplitter

○ MolecularWeightSplitter

○ MaxMinSplitter

○ ButinaSplitter

○ FingerprintSplitter

Let’s take a look how different splitters work.

RandomSplitter

This is one of the simplest splitters. It just selects samples for the training, validation, and test sets in a completely random way.

Didn’t we just say that’s a bad idea? Well, it depends on your data. If every sample is truly independent of every other, then this is just as good a way as any to split the data. There is no universally best choice of splitter. It all depends on your particular dataset, and for some datasets this is a fine choice.

RandomStratifiedSplitter

Some datasets are very unbalanced: only a tiny fraction of all samples are positive. In that case, random splitting may sometimes lead to the validation or test set having few or even no positive samples for some tasks. That makes it unable to evaluate performance.

RandomStratifiedSplitter addresses this by dividing up the positive and negative samples evenly. If you ask for a 80/10/10 split, the validation and test sets will contain not just 10% of samples, but also 10% of the positive samples for each task.

ScaffoldSplitter

This splitter tries to address the problem discussed above where many molecules are very similar to each other. It identifies the scaffold that forms the core of each molecule, and ensures that all molecules with the same scaffold are put into the same dataset. This is still not a perfect solution, since two molecules may have different scaffolds but be very similar in other ways, but it usually is a large improvement over random splitting.

ButinaSplitter

This is another splitter that tries to address the problem of similar molecules. It clusters them based on their molecular fingerprints, so that ones with similar fingerprints will tend to be in the same dataset. The time required by this splitting algorithm scales as the square of the number of molecules, so it is mainly useful for small to medium sized datasets.

SpecifiedSplitter

This splitter leaves everything up to the user. You tell it exactly which samples to put in each dataset. This is useful when you know in advance that a particular splitting is appropriate for your data.

An example is temporal splitting. Consider a research project where you are continually generating and testing new molecules. As you gain more data, you periodically retrain your model on the steadily growing dataset, then use it to predict results for other not yet tested molecules. A good way of validating whether this works is to pick a particular cutoff date, train the model on all data you had at that time, and see how well it predicts other data that was generated later.

TaskSplitter

Provides a simple interface for splitting datasets task-wise.

For some learning problems, the training and test datasets should have different tasks entirely. This is a different paradigm from the usual Splitter, which ensures that split datasets have different data points, not different tasks. This method improves multi-task learning and problem decomposition situations by enhancing their efficiency and performance.

SingletaskStratifiedSplitter

Another way of splitting data, particularly for classification tasks with imbalanced class distributions is the single-task stratified splitter. The single-task stratified splitter maintains the class distribution in the original dataset across training, validation and test sets. This is crucial when working with imbalanced datasets where some classes may be underrepresented.

FingerprintSplitter

Class for doing data splits based on the Tanimoto similarity(Tanimoto similarity measures overlap between two sets succinctly) between ECFP4 fingerprints(ECFP4 fingerprints encode unique parts of molecules for efficient comparison).

This class tries to split the data such that the molecules in each dataset are as different as possible from the ones in the other datasets. This makes it a very stringent test of models. Predicting the test and validation sets may require extrapolating far outside the training data.It splits molecular datasets using Tanimoto similarity scores calculated from ECFP4 fingerprints. ECFP4, based on Morgan fingerprints, encodes molecular substructures.

MolecularWeightSplitter

Another splitter performs data splits based on molecular weight

Effect of Using Different Splitters

Let’s look at an example. We will load the Tox21 toxicity dataset using random, fingerprint, scaffold, and Butina splitting. For each one we train a model and evaluate it on the training and test sets.

1
2
3
4
5
6
7
8
9
10
11
12
13
14

import deepchem as dc

splitters = ['random', 'scaffold', 'butina', 'fingerprint']
metric = dc.metrics.Metric(dc.metrics.roc_auc_score)

for splitter in splitters:
    tasks, datasets, transformers = dc.molnet.load_tox21(featurizer='ECFP', splitter=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()

All of them produce very similar performance on the training set, but the random splitter has much higher performance on the test set. Scaffold splitting has a lower test set score, and Butina splitting is even lower. Does that mean random splitting is better? No! It means random splitting doesn’t give you an accurate measure of how well your model works. Because the test set contains lots of molecules that are very similar to ones in the training set, it isn’t truly independent. It makes the model appear to work better than it really does. Scaffold splitting and Butina splitting give a bette indication of what you can expect on independent data in the future.

Congratulations! Time to join the Community!

Congratulations on completing this tutorial notebook! If you enjoyed working through the tutorial, and want to continue working with DeepChem, we encourage you to finish the rest of the tutorials in this series. You can also help the DeepChem community in the following ways:

Star DeepChem on GitHub

This helps build awareness of the DeepChem project and the tools for open source drug discovery that we’re trying to build.

Join the DeepChem Discord

The DeepChem discord hosts a number of scientists, developers, and enthusiasts interested in deep learning for the life sciences. Join the conversation!

Citing This Tutorial

If you found this tutorial useful please consider citing it using the provided BibTeX.

1
2
3
4
5

In [5]: @manual{Intro8,
title={Working With Splitters},
organization={DeepChem},
author={Eastman, Peter, Mohapatra, Bibhusundar and Ramsundar, Bharath},
howpublished = {\url{https://github.com/deepchem/deepchem/blob/master/examples/tutorials/Working_With_Splitter:year={2021},}

Advanced Model Training

In the tutorials so far we have followed a simple procedure for training models: load a dataset, create a model, call fit() , evaluate it, and call ourselves done. That’s fine for an example, but in real machine learning projects the process is usually more complicated. In this tutorial we will look at a more realistic workflow for training a model.

Colab

This tutorial and the rest in this sequence can be done in Google colab. If you’d like to open this notebook in colab, you can use the following link.

Open in Colab

Setup

To run DeepChem within Colab, you’ll need to run the following installation commands. You can of course run this tutorial locally if you prefer. In that case, don’t run these cells since they will download and install DeepChem in your local machine again.

1

In [ ]: !pip install --pre deepchem import deepchem deepchem.__version__

Hyperparameter Optimization

Let’s start by loading the HIV dataset. It classifies over 40,000 molecules based on whether they inhibit HIV replication.

1
2

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

Now let’s train a model on it. We will use a MultitaskClassifier , which is just a stack of dense layers. But that still leaves a lot of options. How many layers should there be, and how wide should each one be? What dropout rate should we use? What learning rate?

These are called hyperparameters. The standard way to select them is to try lots of values, train each model on the training set, and evaluate it on the validation set. This lets us see which ones work best.

You could do that by hand, but usually it’s easier to let the computer do it for you. DeepChem provides a selection of hyperparameter optimization algorithms, which are found in the dc.hyper package. For this example we’ll use GridHyperparamOpt , which is the most basic method. We just give it a list of options for each hyperparameter and it exhaustively tries all combinations of them.

The lists of options are defined by a dict that we provide. For each of the model’s arguments, we provide a list of values to try. In this example we consider three possible sets of hidden layers: a single layer of width 500, a single layer of width 1000, or two layers each of width 1000. We also consider two dropout rates (20% and 50%) and two learning rates (0.001 and 0.0001).

1
2
3
4
5
6
7
8
9
10

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() returns three arguments: the best model it found, the hyperparameters for that model, and a full listing of the validation score for every model. Let’s take a look at the last one.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

In [4]: 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)

Step 1000 validation: roc_auc_score=0.759757
Step 2000 validation: roc_auc_score=0.770685
Step 3000 validation: roc_auc_score=0.771588
Step 4000 validation: roc_auc_score=0.777862
Step 5000 validation: roc_auc_score=0.773894
Step 6000 validation: roc_auc_score=0.763762
Step 7000 validation: roc_auc_score=0.766361
Step 8000 validation: roc_auc_score=0.767026
Step 9000 validation: roc_auc_score=0.761239
Step 10000 validation: roc_auc_score=0.761279
Step 11000 validation: roc_auc_score=0.765363
Step 12000 validation: roc_auc_score=0.769481
Step 13000 validation: roc_auc_score=0.768523
Step 14000 validation: roc_auc_score=0.761306
Step 15000 validation: roc_auc_score=0.77397
Step 16000 validation: roc_auc_score=0.764848

Out[4]: 0.8040038299560547

1
2
3
4
5

{'_dropouts_0.200000_layer_sizes[500]_learning_rate_0.001000_n_features_1024_n_tasks_1': 0.759624393738977, '_dropouts_0.200000_layer_sizes[500]_learning_rate_0.000100_n_features_1024_n_tasks_1': 0.7680791323731138, '_dropouts_0.500000_layer_sizes[500]_learning_rate_0.001000_n_features_1024_n_tasks_1': 0.7623870149911817, '_dropouts_0.500000_layer_sizes[500]_learning_rate_0.000100_n_features_1024_n_tasks_1': 0.7552282358416618, '_dropouts_0.200000_layer_sizes[1000]_learning_rate_0.001000_n_features_1024_n_tasks_1': 0.7689915858318636, '_dropouts_0.200000_layer_sizes[1000]_learning_rate_0.000100_n_features_1024_n_tasks_1': 0.7619292572996277, '_dropouts_0.500000_layer_sizes[1000]_learning_rate_0.001000_n_features_1024_n_tasks_1': 0.7641491524593376, '_dropouts_0.500000_layer_sizes[1000]_learning_rate_0.000100_n_features_1024_n_tasks_1': 0.7609877155594749, '_dropouts_0.200000_layer_sizes[1000, 1000]_learning_rate_0.001000_n_features_1024_n_tasks_1': 0.7707169802077
21,
' _dropouts_ 0.200000_layer_sizes[1000, 1000] _learning_rate_ 0.                                                                                                    <nl>
329,
' _dropouts_ 0.500000_layer_sizes[1000, 100O] _learning_rate_  {   } . {   } . {   } . {   } . {   } . {   } . {   } . {   } . {   } . {   } . {   } . {   } . {   } . {   } . {   } . {   } . {   } . {   } . {   } . {   } . {   } . {   } . {   } . {   } . {   } . {   }

We can see a few general patterns. Using two layers with the larger learning rate doesn’t work very well. It seems the deeper model requires a smaller learning rate. We also see that 20% dropout usually works better than 50%. Once we narrow down the list of models based on these observations, all the validation scores are very close to each other, probably close enough that the remaining variation is mainly noise. It doesn’t seem to make much difference which of the remaining hyperparameter sets we use, so let’s arbitrarily pick a single layer of width 1000 and learning rate of 0.0001.

Early Stopping

There is one other important hyperparameter we haven’t considered yet: how long we train the model for. GridHyperparamOpt trains each for a fixed, fairly small number of epochs. That isn’t necessarily the best number.

You might expect that the longer you train, the better your model will get, but that isn’t usually true. If you train too long, the model will usually start overfitting to irrelevant details of the training set. You can tell when this happens because the validation set score stops increasing and may even decrease, while the score on the training set continues to improve.

Fortunately, we don’t need to train lots of different models for different numbers of steps to identify the optimal number. We just train it once, monitor the validation score, and keep whichever parameters maximize it. This is called “early stopping”. DeepChem’s ValidationCallback class can do this for us automatically. In the example below, we have it compute the validation set’s ROC AUC every 1000 training steps. If you add the save_dir argument, it will also save a copy of the best model parameters to disk.

Learning Rate Schedules

In the examples above we use a fixed learning rate throughout training. In some cases it works better to vary the learning rate during training. To do this in DeepChem, we simply specify a LearningRateSchedule object instead of a number for the learning_rate argument. In the following example we use a learning rate that decreases exponentially. It starts at 0.0002, then gets multiplied by 0.9 after every 1000 steps.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

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)

Step 1000 validation: roc_auc_score=0.736547
Step 2000 validation: roc_auc_score=0.758979
Step 3000 validation: roc_auc_score=0.768361
Step 4000 validation: roc_auc_score=0.764898
Step 5000 validation: roc_auc_score=0.775253
Step 6000 validation: roc_auc_score=0.779898
Step 7000 validation: roc_auc_score=0.76991
Step 8000 validation: roc_auc_score=0.771515
Step 9000 validation: roc_auc_score=0.773796
Step 10000 validation: roc_auc_score=0.776977
Step 11000 validation: roc_auc_score=0.778866
Step 12000 validation: roc_auc_score=0.777066
Step 13000 validation: roc_auc_score=0.77616
Step 14000 validation: roc_auc_score=0.775646
Step 15000 validation: roc_auc_score=0.772785
Step 16000 validation: roc_auc_score=0.769975

Out[5]: 0.22854619979858398

Congratulations! Time to join the Community!

Congratulations on completing this tutorial notebook! If you enjoyed working through the tutorial, and want to continue working with DeepChem, we encourage you to finish the rest of the tutorials in this series. You can also help the DeepChem community in the following ways:

Star DeepChem on GitHub

This helps build awareness of the DeepChem project and the tools for open source drug discovery that we’re trying to build.

Join the DeepChem Discord

The DeepChem Discord hosts a number of scientists, developers, and enthusiasts interested in deep learning for the life sciences. Join the conversation!

Citing This Tutorial

If you found this tutorial useful please consider citing it using the provided BibTeX.

1
2
3
4
5
6
7

@manual{Intro9,
title={Advanced Model Training},
organization={DeepChem},
author={Eastman, Peter and Ramsundar, Bharath},
howpublished = {\url{https://github.com/deepchem/deepchem/blob/master/examples/tutorials/Advanced_Model_Training},
year={2021},
}

Creating a High Fidelity Dataset from Experimental Data

In this tutorial, we will look at what is involved in creating a new Dataset from experimental data. As we will see, the mechanics of creating the Dataset object is only a small part of the process. Most real datasets need significant cleanup and QA before they are suitable for training models.

Colab

This tutorial and the rest in this sequence are designed to be done in Google colab. If you’d like to open this notebook in colab, you can use the following link.

Open in Colab

!pip install –pre deepchemimport deepchemdeepchem.__version

Working With Data Files

Suppose you were given data collected by an experimental collaborator. You would like to use this data to construct a machine learning model.

How do you transform this data into a dataset capable of creating a useful model?

Building models from novel data can present several challenges. Perhaps the data was not recorded in a convenient manner. Additionally, perhaps the data contains noise. This is a common occurrence with, for example, biological assays due to the large number of external variables and the difficulty and cost associated with collecting multiple samples. This is a problem because you do not want your model to fit to this noise.

Hence, there are two primary challenges:

Parsing data

De-noising data

In this tutorial, we will walk through an example of curating a dataset from an excel spreadsheet of experimental drug measurements. Before we dive into this example though, let’s do a brief review of DeepChem’s input file handling and featurization capabilities.

Input Formats

DeepChem supports a whole range of input files. For example, accepted input formats include .csv, .sdf, .fasta, .png, .tif and other file formats. The loading for a particular file format is governed by the Loader class associated with that format. For example, to load a .csv file we use the CSVLoader class. Here’s an example of a .csv file that fits the requirements of CSVLoader .

1. A column containing SMILES strings.

2. A column containing an experimental measurement.

3. (Optional) A column containing a unique compound identifier.

Here’s an example of a potential input file.

Compound ID measured log solubility in mols per litre smiles

Here the “smiles” column contains the SMILES string, the “measured log solubility in mols per litre” contains the experimental measurement, and “Compound ID” contains the unique compound identifier.

Data Featurization

Most machine learning algorithms require that input data form vectors. However, input data for drug-discovery datasets routinely come in the form of lists of molecules and associated experimental readouts. To load the data, we use a subclass of dc.data.DataLoader such as dc.data.CSVLoader or dc.data.SDFLoader . Users can subclass dc.data.DataLoader to load arbitrary file formats. All loaders must be passed a dc.feat.Featurizer object, which specifies how to transform molecules into vectors. DeepChem provides a number of different subclasses of dc.feat.Featurizer .

Parsing data

In order to read in the data, we will use the pandas data analysis library.

In order to convert the drug names into smiles strings, we will use pubchempy. This isn’t a standard DeepChem dependency, but you can install this library with conda install pubchempy

1
2
3
4

In [ ]: !conda install pubchempy
In [1]: import os
import pandas as pd
from pubchempy import get_cids, get_compounds

Pandas is magic but it doesn’t automatically know where to find your data of interest. You likely will have to look at it first using a GUI.

We will now look at a screenshot of this dataset as rendered by LibreOffice.

To do this, we will import Image and os.

1
2
3
4
5

In [2]: 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))

We see the data of interest is on the second sheet, and contained in columns “TA ID”, “N #1 (%)”, and “N #2 (%)”.

Additionally, it appears much of this spreadsheet was formatted for human readability (multicolumn headers, column labels with spaces and symbols, etc.). This makes the creation of a neat dataframe object harder. For this reason we wil cut everything that is unnecesary or inconvenient.

Note that the actual row headers are stored in row 1 and not 0 above.

1
2
3
4
5
6
7
8
9

In [6]: # 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']

1

In [7]: # preview cleaner dataframe raw_data.loc[raw_data.index[:5]]

This formatting is closer to what we need.

Now, let’s take the drug names and get smiles strings for them (format needed for DeepChem).

1

In [8]: drugs = raw_data['drug'].values

For many of these, we can retreive the smiles string via the canonical_smiles attribute of the get_compounds object (using pubchempy )

1

In [9]: get_compounds(drugs[1], 'name')

[Compound(5281078)]

1

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

1

Out[10]: 'CC1=C2C0C(=0)C2=C(C(=C10C)CC=C(C)CCC(=0)OCCN3CC0CC3)0'

However, some of these drug names have variables spaces and symbols (·, (±), etc.), and names that may not be readable by pubchempy.

For this task, we will do a bit of hacking via regular expressions. Also, we notice that all ions are written in a shortened form that will need to be expanded. For this reason we use a dictionary, mapping the shortened ion names to versions recognizable to pubchempy.

Unfortunately you may have several corner cases that will require more hacking.

1
2
3
4
5
6
7
8
9
10

In [11]: import re
ion_replacements = {
    'HBr': 'hydrobromide',
    '2Br': 'dibromide',
    'Br': 'bromide',
    'HCl': 'hydrochloride',
    '2H2O': 'dihydrate',
    'H2O': 'hydrate',
    'Na': 'sodium'
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

}
ion_keys = ['H20', 'HBr', 'HCl', '2Br', '2H20', '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

Now let’s actually convert all these compounds to smiles. This conversion will take a few minutes so might not be a bad spot to go grab a coffee or tea and take a break while this is running! Note that this conversion will sometimes fail so we’ve added some error handling to catch these cases below.

1
2
3
4
5
6
7
8
9

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

1
2
3

In [13]: 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)

1
2

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

1
2
3
4
5
6

label drug n1 n2
0 2 CC1(C(N2C(S1)C(C2=O)NC(=O)COC3=CC=CC=C3)C(=O)[... -10.404 -18.1929
1 3 CC1=C2COC(=O)C2=C(C(=C1OC)CC=C(C)CCC(=O)OCCN3C... -12.4453 -11.7175
2 4 CC1=CC(=CC(=C1)OCC2CNC(=O)O2)C -8.65572 -17.7753
3 5 COC1=C(C=C2C(=C1)C(=NC(=N2)N3CCN(CC3)C(=O)C4CC... -11.5048 16.0825
4 6 CC(C)N1C2=CC=CC=C2C(=C1C=CC(CC(CC(=O)[O-])O)O)... -11.1354 -14.553

Hooray, we have mapped each drug name to its corresponding smiles code.

Now, we need to look at the data and remove as much noise as possible.

De-noising data

In machine learning, we know that there is no free lunch. You will need to spend time analyzing and understanding your data in order to frame your problem and determine the appropriate model framework. Treatment of your data will depend on the conclusions you gather from this process.

Questions to ask yourself:

What are you trying to accomplish?

What is your assay?

What is the structure of the data?

Does the data make sense?

What has been tried previously?

For this project (respectively):

I would like to build a model capable of predicting the affinity of an arbitrary small molecule drug to a particular ion channel protein

For an input drug, data describing channel inhibition

A few hundred drugs, with n=2

Will need to look more closely at the dataset*

Nothing on this particular protein

*This will involve plotting, so we will import matplotlib and seaborn. We will also need to look at molecular structures, so we will import rdkit. We will also use the seaborn library which you can install with conda install seaborn

1
2
3
4
5
6
7
8
9
10
11
12

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

Our goal is to build a small molecule model, so let’s make sure our molecules are all small. This can be approximated by the length of each smiles string.

1
2
3
4
5

In [16]: 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')

Some of these look rather large, len(smiles) > 150. Let’s see what they look like.

1
2
3
4
5
6
7
8

In [17]: # 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=6)

As suspected, these are not small molecules, so we will remove them from the dataset. The argument here is that these molecules could register as inhibitors simply because they are large. They are more likely to sterically block the channel, rather than diffuse inside and bind (which is what we are interested in).

The lesson here is to remove data that does not fit your use case.

1
2
3
4

In [18]: # drop large molecules
smiles_data = smiles_data[~smiles_data['drug'].isin(long_smiles)]
Now, let's look at the numerical structure of the dataset.
First, check for NaNs.

1

Out[21]: Text(0, 0.5, 'probability')

1

Out[20]: <seaborn.axisgrid.JointGrid at 0x14c4e37d0>

1
2

In [19]: nan_rows = smiles_data[smiles_data.isnull().T.any().T]
nan_rows[['n1', 'n2']]

I don’t trust n=1, so I will throw these out.

Then, let’s examine the distribution of n1 and n2.

1
2
3

In [20]: 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)

We see that most of the data is contained in the gaussian-ish blob centered a bit below zero. We see that there are a few clearly active datapoints located in the bottom left, and one on the top right. These are all distinguished from the majority of the data. How do we handle the data in the blob?

Because n1 and n2 represent the same measurement, ideally they would be of the same value. This plot should be tightly aligned to the diagonal, and the pearson correlation coefficient should be 1. We see this is not the case. This helps gives us an idea of the error of our assay.

Let’s look at the error more closely, plotting in the distribution of (n1-n2).

1
2
3
4

In [21]: diff_df = df['n1'] - df['n2']
sns.histplot(diff_df)
plt.xlabel('difference in n')
plt.ylabel('probability')

This looks pretty gaussian, let’s get the 95% confidence interval by fitting a gaussian via scipy, and taking 2*the standard deviation

1
2
3
4

In [22]: from scipy import stats
mean, std = stats.norm.fit(np.asarray(diff_df, dtype=np.float32))
ci_95 = std*2
ci_95

1

Out[22]: 17.75387954711914

Now, I don’t trust the data outside of the confidence interval, and will therefore drop these datapoints from df.

For example, in the plot above, at least one datapoint has n1-n2 > 60. This is disconcerting.

1
2
3

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

Now that data looks much better!

So, let’s average n1 and n2, and take the error bar to be ci_95.

1
2
3
4

In [24]: 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)

Now, let’s look at the sorted data with error bars.

1
2
3

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

1

Out[25]: Text(0, 0.5, 'activity')

Now, let’s identify our active compounds.

In my case, this required domain knowledge. Having worked in this area, and having consulted with professors specializing on this channel, I am interested in compounds where the absolute value of the activity is greater than 25. This relates to the desired drug potency we would like to model.

If you are not certain how to draw the line between active and inactive, this cutoff could potentially be treated as a hyperparameter.

1
2

In [26]: actives = avg_df[abs(avg_df['n']) - ci_95 > 25]['n']
plt.errorbar(np.arange(actives.shape[0]), actives, yerr=ci_95, fmt='o')

1
2
3

In [27]: # summary
print (raw_data.shape, avg_df.shape, len(actives.index))
(430, 5) (392, 3) 6

In summary, we have:

Removed data that did not address the question we hope to answer (small molecules only)

Dropped NaNs

Determined the noise of our measurements

Removed exceptionally noisy datapoints

Identified actives (using domain knowledge to determine a threshold)

Determine model type, final form of dataset, and sanity load

Now, what model framework should we use?

Given that we have 392 datapoints and 6 actives, this data will be used to build a low data one-shot classifier (10.1021/acscentsci.6b00367). If there were datasets of similar character, transfer learning could potentially be used, but this is not the case at the moment.

Let’s apply logic to our dataframe in order to cast it into a binary format, suitable for classification.

1
2
3

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

1

In [29]: avg_df.to_csv('modulators.csv', index=False)

Now, we will convert this dataframe to a DeepChem dataset.

1
2
3
4
5
6

In [30]: 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)

Lastly, it is often advantageous to numerically transform the data in some way. For example, sometimes it is useful to normalize the data, or to zero the mean. This depends in the task at hand.

Built into DeepChem are many useful transformers, located in the deepchem.transformers.transformers base class.

Because this is a classification model, and the number of actives is low, I will apply a balancing transformer. I treated this transformer as a hyperparameter when I began training models. It proved to unambiguously improve model performance.

1
2

In [31]: transformer = dc.trans.BalancingTransformer(dataset=dataset)
    dataset = transformer.transform(dataset)

Now let’s save the balanced dataset object to disk, and then reload it as a sanity check.

1
2

In [32]: dc.utils.save_to_disk(dataset, 'balanced_dataset.joblib')
balanced_dataset = dc.utils.load_from_disk('balanced_dataset.joblib')

Tutorial written by Keri McKiernan (github.com/kmckiern) on September 8, 2016

Congratulations! Time to join the Community!

Congratulations on completing this tutorial notebook! If you enjoyed working through the tutorial, and want to continue working with DeepChem, we encourage you to finish the rest of the tutorials in this series. You can also help the DeepChem community in the following ways:

Star DeepChem on GitHub

This helps build awareness of the DeepChem project and the tools for open source drug discovery that we’re trying to build.

Join the DeepChem Discord

The DeepChem Discord hosts a number of scientists, developers, and enthusiasts interested in deep learning for the life sciences. Join the conversation!

Bibliography

[2] Anderson, Eric, Gilman D. Veith, and David Weininger. “SMILES, a line notation and computerized interpreter for chemical structures.” US Environmental Protection Agency, Environmental Research Laboratory, 1987.

Citing This Tutorial

If you found this tutorial useful please consider citing it using the provided BibTeX.

1
2
3
4
5
6
7

In [ ]: @manual{Intro10,
title={Creating a high fidelity model from experimental data},
organization={DeepChem},
author={Eastman, Peter and Ramsundar, Bharath},
howpublished = {\url{https://github.com/deepchem/deepchem/tree/master/examples/tutorials}},
year={2021},
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

In [3]: 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)

MUV-466 0.9207684040838259
MUV-548 0.7480655561526062
MUV-600 0.9927995701235895
MUV-644 0.9974207415368082
MUV-652 0.7823481998925309
MUV-689 0.6636843990686011
MUV-692 0.6319093677234462
MUV-712 0.7787838079885365
MUV-713 0.7910711087229088
MUV-733 0.4401307540748701
MUV-737 0.34679383843811573
MUV-810 0.9564571019165323
MUV-832 0.9991044241447251
MUV-846 0.7519881783987103
MUV-852 0.8516747268493642
MUV-858 0.5906591438294824
MUV-859 0.5962954008166774

Putting Multitask Learning to Work

This notebook walks through the creation of multitask models on MUV [1]. The goal is to demonstrate how multitask methods can provide improved performance in situations with little or very unbalanced data.

Colab

This tutorial and the rest in this sequence are designed to be done in Google colab. If you’d like to open this notebook in colab, you can use the following link.

1

In [ ]: !pip install --pre deepchem import deepchem deepchem.__version__

The MUV dataset is a challenging benchmark in molecular design that consists of 17 different “targets” where there are only a few “active” compounds per target. There are 93,087 compounds in total, yet no task has more than 30 active compounds, and many have even less. Training a model with such a small number of positive examples is very challenging. Multitask models address this by training a single model that predicts all the different targets at once. If a feature is useful for predicting one task, it often is useful for predicting several other tasks as well. Each added task makes it easier to learn important features, which improves performance on other tasks [2].

To get started, let’s load the MUV dataset. The MoleculeNet loader function automatically splits it into training, validation, and test sets. Because there are so few positive examples, we use stratified splitting to ensure the test set has enough of them to evaluate.

1
2
3
4
5

In [1]: import deepchem as dc
import numpy as np

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

Now let’s train a model on it. We’ll use a MultitaskClassifier, which is a simple stack of fully connected layers.

1
2
3
4

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

1

Out[2]: 0.0004961589723825455

Let’s see how well it does on the test set. We loop over the 17 tasks and compute the ROC AUC for each one.

Not bad! Recall that random guessing would produce a ROC AUC score of 0.5, and a perfect predictor would score 1.0. Most of the tasks did much better than random guessing, and many of them are above 0.9.

Congratulations! Time to join the Community!

Congratulations on completing this tutorial notebook! If you enjoyed working through the tutorial, and want to continue working with DeepChem, we encourage you to finish the rest of the tutorials in this series. You can also help the DeepChem community in the following ways:

Star DeepChem on GitHub

This helps build awareness of the DeepChem project and the tools for open source drug discovery that we’re trying to build.

Join the DeepChem Gitter

The DeepChem Gitter hosts a number of scientists, developers, and enthusiasts interested in deep learning for the life sciences. Join the conversation!

Bibliography

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

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

Tutorial Part 13: Modeling Protein-Ligand Interactions

By Nathan C. Frey | Twitter and Bharath Ramsundar | Twitter

In this tutorial, we’ll walk you through the use of machine learning and molecular docking methods to predict the binding energy of a protein-ligand complex. Recall that a ligand is some small molecule which interacts (usually noncovalently) with a protein. Molecular docking performs geometric calculations to find a “binding pose” with a small molecule interacting with a protein in a suitable binding pocket (that is, a region on the protein which has a groove in which the small molecule can rest).

The structure of proteins can be determined experimentally with techniques like Cryo-EM or X-ray crystallography. This can be a powerful tool for structure-based drug discovery. For more info on docking, read the AutoDock Vina paper and the deepchem.dock documentation. There are many graphical user and command line interfaces (like AutoDock) for performing molecular docking. Here, we show how docking can be performed programmatically with DeepChem, which enables automation and easy integration with machine learning pipelines.

As you work through the tutorial, you’ll trace an arc including

1. Loading a protein-ligand complex dataset (PDBbind)

2. Performing programmatic molecular docking

3. Featurizing protein-ligand complexes with interaction fingerprints

4. Fitting a random forest model and predicting binding affinities

To start the tutorial, we’ll use a simple pre-processed dataset file that comes in the form of a gzipped file. Each row is a molecular system, and each column represents a different piece of information about that system. For instance, in this example, every row reflects a protein-ligand complex, and the following columns are present: a unique complex identifier; the SMILES string of the ligand; the binding affinity (Ki) of the ligand to the protein in the complex; a Python list of all lines in a PDB file for the protein alone; and a Python list of all lines in a ligand file for the ligand alone.

Colab

This tutorial and the rest in this sequence are designed to be done in Google colab. If you’d like to open this notebook in colab, you can use the following link.

Open in Colab

Setup

To run DeepChem within Colab, you’ll need to run the following cell of installation commands. This will take about 5 minutes to run to completion and install your environment.

1
2
3

!pip install -q condacolab
import condacolab
condacolab.install()

WARNING: Running pip as the ‘root’ user can result in broken permissions and conflicting behaviour with the syst em package manager. It is recommended to use a virtual environment instead: https://pip.pypa.io/warnings/venv Everything looks OK!

1

In [2]: !conda install -c conda-forge openmm

Collecting package metadata (current_repodata.json): - \ | 7 - \ | / - \ | - \ - \ 一 / - \ | - \ | / - / - \ | / - | / - \ | / - \ | / \ | / - \ / - \ | / - \ / - \ / - \ / - / - / - \ / - \ | / \ | / - \ | / - \ | / - \ | / - \ | / - \ / - \ | / - \ / \ - \ / - / - / \ / - \ | / \ | / 1 \ / - \ | - \ | - \ | / - \ | / - \ | / - \ | \ | - \ | - \ / - \ / - \ - \ / - | / - \ / done Solving environment: \ | / - \ | 7 / - \ | - \ / - \ | - \ / - \ | / done

# All requested packages already installed.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/
Requirement already satisfied: deepchem in /usr/local/lib/python3.10/site-packages (2.7.1)
Requirement already satisfied: numpy>=1.21 in /usr/local/lib/python3.10/site-packages (from deepchem) (1.24.3)
Requirement already satisfied: scikit-learn in /usr/local/lib/python3.10/site-packages (from deepchem) (1.2.2)
Requirement already satisfied: joblib in /usr/local/lib/python3.10/site-packages (from deepchem) (1.2.0)
Requirement already satisfied: scipy<1.9 in /usr/local/lib/python3.10/site-packages (from deepchem) (1.8.1)
Requirement already satisfied: pandas in /usr/local/lib/python3.10/site-packages (from deepchem) (2.0.1)
Requirement already satisfied: rdkit in /usr/local/lib/python3.10/site-packages (from deepchem) (2023.3.1)
Requirement already satisfied: python-dateutil>=2.8.2 in /usr/local/lib/python3.10/site-packages (from pandas->deepchem) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.10/site-packages (from pandas->deepchem) (2023.3)
Requirement already satisfied: tzdata>=2022.1 in /usr/local/lib/python3.10/site-packages (from pandas->deepchem) (2023.3)
Requirement already satisfied: Pillow in /usr/local/lib/python3.10/site-packages (from rdkit->deepchem) (9.5.0)
Requirement already satisfied: threadpoolctl>=2.0.0 in /usr/local/lib/python3.10/site-packages (from scikit-learn->deepchem) (3.1.0)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.10/site-packages (from python-dateutil>=2.8.2->pandas->deepchem) (1.16.0)
WARNING: Running pip as the 'root' user can result in broken permissions and conflicting behaviour with the system package manager. It is recommended to use a virtual environment instead: https://pip.pypa.io/warnings/venv

1

In [4]: !conda install -c conda-forge pdbfixer

1
2
3
4
5

Collecting package metadata (current_repodata.json): - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ |
| / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ |
| / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ |
| / - \ | / - \ | / - \ | / - \$ | / - done
Solving environment: | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ done

# All requested packages already installed.

1

In [5]: !conda install -c conda-forge vina

1
2
3
4
5
6

Collecting package metadata (current_repodata.json): - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ |
| / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ |
| / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ | / - \ |
| / - \ | / - \ | / - \ | / - \$ | / - \$ | / - \$ | / - \$ | / - \$ | / - \$ | / - \$ | / - \$ | / - \$ | / - \$ | / - \$ |
| / - \ | / - \ | / - \ | / - \$ | / - \$ | / - \$ done
Solving environment: / - \ | / - \ | / - \ | / - \$ | / - \$ | / - \$ | / - \$ | / - \$ | / - \$ | / - \$ | / - \$ | / - \$ | / done

1

# All requested packages already installed.

Protein-ligand complex data

It is really helpful to visualize proteins and ligands when doing docking. Unfortunately, Google Colab doesn’t currently support the Jupyter widgets we need to do that visualization. Install MDTraj and nglview on your local machine to view the protein-ligand complexes we’re working with.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

In [6]: !pip install -q mdtraj nglview
# !jupyter-nbextension enable nglview --py --sys-prefix # for jupyter notebook
# !jupyter labextension install nglview-js- widgets # for jupyter lab

WARNING: Running pip as the 'root' user can result in broken permissions and conflicting behaviour with the system package manager. It is recommended to use a virtual environment instead: https://pip.pypa.io/warnings/venv

In [8]: 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

Skipped loading modules with pytorch-geometric dependency, missing a dependency. No module named 'torch_geometric'
Skipped loading modules with pytorch-geometric dependency, missing a dependency. cannot import name 'DMPNN' from 'deepchem.models.torch_models' (/usr/local/lib/python3.10/site-packages/deepchem/models/torch_models/_init_.py)
Skipped loading modules with pytorch-lightning dependency, missing a dependency. No module named 'pytorch_lightning'
Skipped loading some Jax models, missing a dependency. No module named 'haiku'

To illustrate the docking procedure, here we’ll use a csv that contains SMILES strings of ligands as well as PDB files for the ligand and protein targets from PDBbind. Later, we’ll use the labels to train a model to predict binding affinities. We’ll also show how to download and featurize PDBbind to train a model from scratch.

1
2
3
4
5

<timed exec>:7: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.
Warning: importing 'simtk.openmm' is deprecated. Import 'openmm' instead.
3cyx 1510
CPU times: user 2.04 s, sys: 157 ms, total: 2.2 s
Wall time: 4.32 s

1
2
3
4
5
6
7
8
9
10
11

In [9]:
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 sponsor')

raw_dataset = load_from_disk(dataset_file)
raw_dataset = raw_dataset[['pdb_id', 'smiles', 'label']]

1

In [10]: raw_dataset.head(2)

1

Out[10]: pdb_id smiles label 0 2d3u CC1CCCCC1S(O)(O)NC1CC(C2CCC(CN)CC2)SC1C(O)O 6.92 1 3cyx CC(C)(C)NC(O)C1CC2CCCCCC2C[NH+]1CC(O)C(CC1CCCCCC... 8.00

Fixing PDB files

Next, let’s get some PDB protein files for visualization and docking. We’ll use the PDB IDs from our raw_dataset and download the pdb files directly from the Protein Data Bank using pdbfixer . We’ll also sanitize the structures with RDKit. This ensures that any problems with the protein and ligand files (non-standard residues, chemical validity, etc.) are corrected. Feel free to modify these cells and pdbids to consider new protein-ligand complexes. We note here that PDB files are complex and human judgement is required to prepare protein structures for docking. DeepChem includes a number of docking utilites to assist you with preparing protein files, but results should be inspected before docking is attempted.

1
2
3

In [11]: from openmm.app import PDBFile
from pdbfixer import PDBFixer
from deepchem.utils.vina_utils import prepare_inputs

1
2
3
4
5
6
7
8
9
10
11
12
13

In [13]: %%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))

Visualization

If you’re outside of Colab, you can expand these cells and use MDTraj and nglview to visualize proteins and ligands.

1
2
3

In [14]: import mdtraj as md
import nglview
from IPython.display import display, Image

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

1
2

In [15]: protein_mdtraj = md.load_pdb('3cyx.pdb')
ligand_mdtraj = md.load_pdb('ligand_3cyx.pdb')

We’ll use the convenience function nglview.show_mdtraj in order to view our proteins and ligands. Note that this wil only work if you uncommented the above cell, installed nglview, and enabled the necessary notebook extensions.

v = nglview.show_mdtraj(ligand_mdtraj)

display(v) # interactive view outside Colab NGLWidget()

Now that we have an idea of what the ligand looks like, let’s take a look at our protein:

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

NGLWidget()

Molecular Docking

Ok, now that we’ve got our data and basic visualization tools up and running, let’s see if we can use molecular docking to estimate the binding affinities between our protein ligand systems.

There are three steps to setting up a docking job, and you should experiment with different settings. The three things we need to specify are 1) how to identify binding pockets in the target protein; 2) how to generate poses (geometric configurations) of a ligand in a binding pocket; and 3) how to “score” a pose. Remember, our goal is to identify candidate ligands that strongly interact with a target protein, which is reflected by the score.

DeepChem has a simple built-in method for identifying binding pockets in proteins. It is based on the convex hull method. The method works by creating a 3D polyhedron (convex hull) around a protein structure and identifying the surface atoms of the protein as the ones closest to the convex hull. Some biochemical properties are considered, so the method is not purely geometrical. It has the advantage of having a low computational cost and is good enough for ou purposes.

the hood, allowing us to get up and running generating poses quickly.

1

In [20]: vpg = dc.dock.pose_generation.VinaPoseGenerator()

We could specify a pose scoring function from deepchem.dock.pose_scoring , which includes things like repulsive and hydrophobic interactions and hydrogen bonding. Vina will take care of this, so instead we’ll allow Vina to compute scores for poses.

1
2
3
4
5
6
7
8
9
10
11
12
13
14

In [7]: !mkdir -p vina_test

In [21]: %%time
complexes, scores = vpg.generate_poses(molecular_complex=('3cyx.pdb', 'ligand_3cyx.pdb'), # protein-ligand file
    out_dir='vina_test',
    generate_scores=True
)

CPU times: user 41min 4s, sys: 21.9 s, total: 41min 26s
Wall time: 28min 32s

/usr/local/lib/python3.10/site-packages/vina/vina.py:260: DeprecationWarning: `np.int` is a deprecated alias for the builtin `int`. To silence this warning, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
self._voxels = np.ceil(np.array(box_size) / self._spacing).astype(np.int)

We used the default value for num_modes when generating poses, so Vina will return the 9 lowest energy poses it found in units of kcal/mol .

1

In [22]: scores

1

Out[22]: [-9.484, -9.405, -9.195, -9.151, -8.9, -8.696, -8.687, -8.633, -8.557]

Can we view the complex with both protein and ligand? Yes, but we’ll need to combine the molecules into a single RDkit molecule.

1

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

Let’s now visualize our complex. We can see that the ligand slots into a pocket of the protein.

1
2

In [24]: v = nglview.show_rdkit(complex_mol) display(v)
NGLWidget()

Now that we understand each piece of the process, we can put it all together using DeepChem’s Docker class. Docker creates a generator that yields tuples of posed complexes and docking scores.

1
2

In [25]: docker = dc.dock.dockering.Docker(pose_generator=vpg)
posed_complex, score = next(docker.dock(molecular_complex=('3cyx.pdb', 'ligand_3cyx.pdb'), use_pose_generator_scores=True))

Modeling Binding Affinity

Docking is a useful, albeit coarse-grained tool for predicting protein-ligand binding affinities. However, it takes some time, especially for large-scale virtual screenings where we might be considering different protein targets and thousands of potential ligands. We might naturally ask then, can we train a machine learning model to predict docking scores? Let’s try and find out!

We’ll show how to download the PDBbind dataset. We can use the loader in MoleculeNet to get the 4852 protein-ligand complexes from the “refined” set or the entire “general” set in PDBbind. For simplicity, we’ll stick with the ~100 complexes we’ve already processed to train our models.

Next, we’ll need a way to transform our protein-ligand complexes into representations which can be used by learning algorithms. Ideally, we’d have neural protein-ligand complex fingerprints, but DeepChem doesn’t yet have a good learned fingerprint of this sort. We do however have well-tuned manual featurizers that can help us with our challenge here.

We’ll make use of two types of fingerprints in the rest of the tutorial, the CircularFingerprint and ContactCircularFingerprint . DeepChem also has voxelizers and grid descriptors that convert a 3D volume containing an arragment of atoms into a fingerprint. These featurizers are really useful for understanding protein-ligand complexes since they allow us to translate complexes into vectors that can be passed into a simple machine learning algorithm. First, we’ll create circular fingerprints. These convert small molecules into a vector of fragments.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

In [26]: pdbids = raw_dataset['pdb_id'].values
ligand_smiles = raw_dataset['smiles'].values

In [29]: %%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))

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.
<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.
<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.
<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.
<timed activ:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.
[15:12:02] UFFTYPER: Warning: hybridization set to SP3 for atom 17
<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.
[15:12:04] UFFTYPER: Warning: hybridization set to SP3 for atom 6
<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.
[15:12:06] UFFTYPER: Warning: hybridization set to SP3 for atom 1
[15:12:06] UFFTYPER: Unrecognized atom type: S_5+4 (21)
<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.
<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.
<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.
<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.
[15:12:23] UFFTYPER: Warning: hybridization set to SP3 for atom 20
<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.
<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.
[15:12:23] UFFTYPER: Warning: hybridization set to SP3 for atom 20

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69

ction in deepchem.utils.docking_utils.

<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.

[15:19:01] UFFTYPER: Warning: hybridization set to SP3 for atom 16

[15:19:01] UFFTYPER: Unrecognized atom type: S_5+4 (20)

<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.

[15:19:02] UFFTYPER: Unrecognized atom type: S_5+4 (6)

<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.

<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.

[15:19:05] UFFTYPER: Unrecognized atom type: S_5+4 (6)

<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.

<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.

1hfs failed sanitization

<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.

<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.

<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.

[15:19:22] UFFTYPER: Warning: hybridization set to SP3 for atom 20

<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.

<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.

[15:19:41] Explicit valence for atom # 1800 C, 5, is greater than permitted

[15:19:41] UFFTYPER: Unrecognized atom type: S_5+4 (11)

<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.

[15:19:42] UFFTYPER: Warning: hybridization set to SP3 for atom 11

<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.

<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.

<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.

[13:19:57] UFFTYPER: Warning: hybridization set to SP3 for atom 9

[15:19:57] UFFTYPER: Warning: hybridization set to SP3 for atom 23

<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.

[15:19:59] UFFTYPER: Warning: hybridization set to SP3 for atom 8

[15:19:59] UFFTYPER: Warning: hybridization set to SP3 for atom 12

[15:19:59] UFFTYPER: Warning: hybridization set to SP3 for atom 34

[15:19:59] UFFTYPER: Warning: hybridization set to SP3 for atom 41

<timed exec>:8: DeprecationWarning: Call to deprecated function prepare_inputs. Please use the corresponding function in deepchem.utils.docking_utils.

CPU times: user 4min 9s, sys: 3.31 s, total: 4min 12s

Wall time: 8min 19s

We’ll do some clean up to make sure we have a valid ligand file for every valid protein. The lines here will compare the PDB IDs between the ligand and protein files and remove any proteins that don’t have corresponding ligands.

1
2
3
4

In [30]: # 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')

1

In [31]: len(proteins), len(ligands)

1

Out[31]: (190, 190)

1
2
3

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

1

In [33]: fp_featurizer = dc.feat.CirclearFingerprint(size=2048)

1

In [34]: features = fp_featurizer.fteaturize([Chem.MolFromPDBFile(l) for l in ligands])

1
2

In [35]: dataset = dc.data.NumpyDataset(X=features, y=labels, ids=pdbids)
train_dataset, test_dataset = dc.splits.RandomSplitter().train_test_split(dataset, seed=42)

The convenience loader dc.molnet.load_pdbbind will take care of downloading and featurizing the pdbbind dataset under the hood for us. This will take quite a bit of time and compute, so the code to do it is commented out. Uncomment it and grab a cup of coffee if you’d like to featurize all of PDBbind’s refined set. Otherwise, you can continue with the small dataset we constructed above.

1
2
3
4

In [36]: # # 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')

Now, we’re ready to do some learning!

To fit a deepchem model, first we instantiate one of the provided (or user-written) model classes. In this case, we have a created a convenience class to wrap around any ML model available in Sci-Kit Learn that can in turn be used to interoperate with deepchem. To instantiate an SklearnModel , you will need (a) task_types, (b) model_params, another dict as illustrated below, and (c) a model_instance defining the type of model you would like to fit, in this case a RandomForestRegressor

1
2
3
4
5

In [38]: 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)

value for the test set indicates that the model isn’t producing meaningful outputs. It turns out that predicting binding affinities is hard. This tutorial isn’t meant to show how to create a state-of-the-art model for predicting binding affinities, but it gives you the tools to generate your own datasets with molecular docking, featurize complexes, and train models.

1
2
3
4
5
6
7
8
9
10
11
12
13

In [39]: # use Pearson correlation so metrics are > 0
metric = dc.metrics.Metric(dc.metrics.pearson_r2_score)

evaluator = Evaluator(model, train_dataset, [])
train_r2score = evaluator.compute_model_performance([metric])
print("RF Train set R^2 %f" % (train_r2score["pearson_r2_score"]))

evaluator = Evaluator(model, test_dataset, [])
test_r2score = evaluator.compute_model_performance([metric])
print("RF Test set R^2 %f" % (test_r2score["pearson_r2_score"]))

RF Train set R^2 0.888697
RF Test set R^2 0.007797

We’re using a very small dataset and an overly simplistic representation, so it’s no surprise that the test set performance is quite bad.

1
2
3
4
5
6
7
8

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

Out[40]: [(6.862549999999994, 7.4),
(6.616400000000008, 6.85),
(4.852004999999995, 3.4),
(6.43060000000001, 6.72),
(8.663229999999999, 11.06)]

1

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

[(5.960549999999999, 4.21), (6.051305714285715, 8.7), (5.799900000000003, 6.39), (6.433881666666665, 4.94), (6.7465399999999995, 9.21)]

The protein-ligand complex view.

In the previous section, we featurized only the ligand. This time, let’s see if we can do something sensible with our protein-ligand fingerprints that make use of our structural information. To start with, we need to re-featurize the dataset but using the contact fingerprint this time.

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

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)

[15:21:40] Explicit valence for atom # 3 C, 5, is greater than permitted Mol [H]OC([H])([H])[C@]1([H])O[C@@]2([H])SC([H])([H])(N([H])C([H])([H])C([H])([H])[H])N([H])[C@@]2([H])[C@@]([H] )(O[H])[C@]1([H])O[H] failed sanitization [15:21:47] Explicit valence for atom # 214 O, 3, is greater than permitted Mol CCC@H[C@@H]1NC(=O)CNC(=O)[C@H](CCSC)NC(=O)[C@H](CCC(N)=O)NC(=O)[C@H](CC2=CC=C(O)C=C2)NC(=O)[C@H](CC2=CN C3=C2C=CC=C3)NC(=O)[C@H](CCCNC(=N)N)NC(=O)[C@H](CC(N)=O)NC(=O)[C@H](CC(N)=O)NC(=O)[C@H](CC2=CC=CC=C2)NC(=O)[C@@H ]2CCCN2C(=O)[C@H](CO)NC(=O)[C@H](CCCCN)NC(=O)[C@H](CCSC)NC(=O)[C@H](C(C)C)NC(=O)[C@H](CC2=CC=CC=C2)NC(=O)[C@@H]2 CCCN2C(=O)CNC(=O)CNC(=O)[C@H](CO)NC(=O)[C@H](CC(=O)O)NC(=O)CNC(=O)[C@H](CCC(=O)O)NC(=O)[C@@H](NC(=O)[C@H](C)NC(= O)[C@H](CC(=O)O)NC(=O)CNC(=O)[C@H](CCCNC(=N)N)NC(=O)[C@H](CCCCN)NC(=O)CNC(=O)[C@H](CCC(=O)O)NC(=O)[C@H](CC(=O)O) NC(=O)[C@@H]2CCCN2C(=O)[C@H](CCCCN)NC(=O)[C@H](CC2=CC=C(O)C=C2)NC(=O)CNC(=O)[C@H](C)NC(=O)[C@@H]2CSSCC@H[C@@H](NC(=O)[C@@H]3CCCN3C(=O)[C@H](CCCNC(=N)N)NC(=O)[C@H](CCC(=O)O)NC(=O)[C@@H](NC(=O)[C@@H](NC(=O)[C@@H]3CCCN 3C(=O)[C@H](CC(C)C)NC(=O)[C@H](CC(N)=O)NC(=O)[C@@H](NC(=O)[C@@H](NC(=O)[C@H](CCC(N)=O)NC(=O)[C@H](CC(C)C)NC(=O) C@@HC(C)C)C(C)C)C(C)C)[C@@H](C)CC)C(C)C)C(C)C)C( =O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)NC@@ HC(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(= O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)N2)CSSC[C@@H](C(=O)NC@@HC(=O)NC@@ HC(=O)NC@@HC(=O)NCC(=O)NC@@HC(=O)NC@@HC(=O)NCC(=O)N[C@@H](CC 2=CC=CC=C2)C(=O)NC@@HC(=O)N[C@H](C(=O)NC@@HC(=O)NC@HN[C@@H](CC2=CC=CC=C 2)C(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)N[C@@H](CC2=CNC3=C2C=CC =C3)C(=O)N[C@H](C(=O)NC@@HC(=O)NC@@HC(=O)N[C@H](C(=O)N[C@H](C(=O)NC@@HC(=O)NC@@ HC(=O)NC@@HC(=O)N2C3C(=O)O32)[C@@H](C)CC)C(C)C)[C@@H](C)CC)C(C)C)[C@@H](C)O)NC(=O)CNC (=O)[C@H](CCC(=O)O)NC(=O)CNC(=O)[C@H](CC2=CNC3=C2C=CC=C3)NC(=O)[C@H](CO)NC(=O)[C@H](C(C)C)NC1=O.CCC@H[C@H]( NC(=O)[C@H](CCC(=O)O)NC(=O)[C@H](C)NC(=O)[C@H](CC(=O)O)NC(=O)[C@H](CO)NC(=O)CNC(=O)[C@]12C(CC(=O)O)C13C(C)[C@H]( N)C(=O)N32)C(=O)NCC(=O)NC@@HC(=O)NC@@HC(=O)N1CCC[C@H]1C(=O)NC@@HC(=O)N[C@@H] (CCC(N)=O)C(=O)N[C@H](C(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC( =O)NC@@HC(=O)NC@@HC(=O)N1CCC[C@H]1C(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)N[C@H]1CSSC[C@@H](C(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)N2CCC[C@H]2C(=O)N2CCC[C@H]2C(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)N C@@HC(=O)NC@@HC(=O)N[C@H](C(=O)NC@@HC(=O)NC@@HC(=O)N[C@@H](CC(= O)O)C(=O)NC@@HC(=O)NC@@HC(=O)N[C@H](C(=O)NC@@HC(=O)N[C@H](C(=O)NCC(=O)NC@@H)C(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)N[C@H](C(=O)N[C@H]2CCCN(C(=N)N)C(C(=O )O)C[C@@H](C(=O)NC@@HC(=O)NC@@HC(=O)N[C@H](C(=O)NC@@HC(=O)NC@@HC(= O)N[C@H](C(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)N[C @H](C(=O)NC@@HC(=O)N[C@H](C(=O)NC@@HC(=O)N3CCC[C@H]3C(=O)NC@@HC(= O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)NC@ @HC(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)NC@@H)C(=O)N[C@H](C(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)NC@@H C(=O)NC@@HC(=O)NC@@HC(=O)N3CCC[C@H]3C(=O)N[C@H](C(=O)NC@@HC(=O)NC@@HC(=O) NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)N[C@H](C(=O)NC@@HC(=O)N3CCC[C@H] 3C(=O)N[C@H](C(=O)N[C@@H](CSSCC@HC(=O)NCC(=O)NC@@HC(=O)NC@@H)C(=O)N3CCC[C@H]3C(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)N[C@@H](CCCC N)C(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)N[C@@H] (CCCCN)C(=O)N[C@H](C(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O) NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)N[C@@]34C(=O)O3C4C3=CC=C(O)C=C3)[C@@H](C)O)C(=O)NC@@H C(=O)N3CCC[C@H]3C(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)NC@HNC @@HC(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O)N[C@@H] (C)C(=O)NCC(=O)NC@@HC(=O)NC@@HC(=O)NCC(=O)NC@@HC(=O)NC@HN[C@H] (C(=O)NCC(=O)NC@@HC(=O)NCC(=O)NC@@HC(=O)NC@@HC(=O)NC@@HC(=O) NC@HCCC(=O)O)[C@@H](C)O)C(C)C)[C@@H](C)O)C(C)C)[C@@H](C)CC)C(C)C)[C@@H](C)CC)[C@@H](C)CC)[C@@H](C)CC)C@@ HCC)[C@@H](C)CC)NC(=O)[C@H](CC3=CC=C(O)C=C3)NC2=O)[C@@H](C)O)[C@@H](C)CC)C(C)C)[C@@H](C)O)NC(=O)[C@H](CC2=CN C=N2)NC(=O)[C@H](C)NC(=O)[C@H](C)NC(=O)[C@H]([C@@H](C)O)NC(=O)[C@H](CC(C)C)NC(=O)[C@H](C(C)C)NC(=O)[C@H](CC2=CNC 3=C2C=CC=C3)NC(=O)[C@H](CCCNC(=N)N)NC(=O)[C@H](CC(=O)O)NC(=O)[C@H](CO)NC(=O)[C@H]([C@@H](C)CC)NC(=O)[C@H](CC(C)C )NC(=O)[C@H](CO)NC(=O)[C@H](C)NC(=O)CNC1=O)C(C)C.CCC@H[C@H](NC(=O)[C@H](CCC(=O)O)NC(=O)[C@H](CCC(=O)O)NC(=O )[C@H](CC1=CC=CC=C1)NC(=O)[C@@H](N)CC(=O)O)C(=O)N1CCC[C@H]1C(=O)NC@@HC(=O)NC@@HC(=O)NC@ HCC1CCC(OS(O)(O)O)CC1.CO[C@H]1CC[C@H](S(O)(O)NC@@HC(O)N2CCC[C@H]2C@H)NCC2CCC(C(N)N)CC2)CC1Cl.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O. O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O. O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O.O. O.O.[NaH].[NaH] failed sanitization