In this tutorial I will walk through how to efficiently screen a large compound library with DeepChem (ZINC). Screening a large compound library using machine learning is a CPU bound pleasingly parrellel problem. The actual code examples I will use assume the resources available are a single very big machine (like an AWS c5.18xlarge), but should be readily swappable for other systmes (like a super computing cluster). At a high level what we will do is...
1. Create a Machine Learning Model Over Labeled Data
2. Transform ZINC into "Work-Units"
3. Create an inference script which runs predictions over a "Work-Unit"
4. Load "Work-Unit" into a "Work Queue"
5. Consume work units from "Work Queue"
6. Gather Results
# 1. Train Model On Labelled Data
We are just going to knock out a simple model here. In a real world problem you will probably try several models and do a little hyper parameter searching.
/home/leswing/miniconda3/envs/deepchem/lib/python3.6/site-packages/tensorflow/python/ops/gradients_impl.py:100: UserWarning: Converting sparse IndexedSlices to a dense Tensor of unknown shape. This may consume a large amount of memory.
"Converting sparse IndexedSlices to a dense Tensor of unknown shape. "
/home/leswing/miniconda3/envs/deepchem/lib/python3.6/site-packages/tensorflow/python/ops/gradients_impl.py:100: UserWarning: Converting sparse IndexedSlices to a dense Tensor of unknown shape. This may consume a large amount of memory.
"Converting sparse IndexedSlices to a dense Tensor of unknown shape. "
%% Cell type:markdown id: tags:
# 2. Create Work-Units
1. Download All of ZINC15.
Go to http://zinc15.docking.org/tranches/home and download all non-empty tranches in .smi format.
I found it easiest to download the wget script and then run the wget script.
For the rest of this tutorial I will assume zinc was downloaded to /tmp/zinc.
The way zinc downloads the data isn't great for inference. We want "Work-Units" which a single CPU can execute that takes a resonable amount of time (10 minutes to an hour). To accomplish this we are going to split the zinc data into files each with 500 thousand lines.
4. splits into multiple files in /tmp/zinc/screen that are 500k molecules long
%% Cell type:markdown id: tags:
## 3. Creat Inference Script
Now that we have work unit we need to construct a program which ingests a work unit and logs the result. It is important that the logging mechanism is thread safe!
For this example we will get the work unit via a file-path, and log the result to a file.
An easy extensions to distribute over multiple computers would be to get the work unit via a url, and log the results to a distributed queue.
We are going to use a flat file as our distribution mechanism. It will be a bash script calling our inference script for every work unit. If you are at an academic institution this would be queing your jobs in pbs/qsub/slurm. An option for cloud computing would be rabbitmq or kafka.
%% Cell type:code id: tags:
``` python
importos
work_units=os.listdir('/tmp/zinc/screen')
withopen('/tmp/zinc/work_queue.sh','w')asfout:
fout.write("#!/bin/bash\n")
forwork_unitinwork_units:
full_path=os.path.join('/tmp/zinc',work_unit)
fout.write("python inference.py %s"%full_path)
```
%% Cell type:markdown id: tags:
# 5. Consume work units from "distribution mechanism"
We will consume work units from our work queue using a very simple Process Pool. It takes lines from our "Work Queue" and runs them, running as many processes in parrallel as we have cpus. If you are using a supercomputing cluster system like pbs/qsub/slurm it will take care of this for you. The key is to use one CPU per work unit to get highest throughput. We accomplish that here using the linux utility "taskset".
Using an c5.18xlarge on aws this will finish overnight.
Since we logged our results to \*_out.smi we now need to gather all of them up and sort them by our predictions. The resulting file wile be > 40GB. To analyze the data further you can use dask, or put the data in a rdkit postgres cartridge.
Here I show how to join the and sort the data to get the "best" results.
The screen seems to favor molecules with one or multiple sulfur trioxides. The top scoring molecules also have low diversity. When creating a "buy list" we want to optimize for more things than just activity, for instance diversity and drug like MPO.
In this example, we will use a `SeqToSeq` model to generate fingerprints for classifying molecules. This is based on the following paper, although some of the implementation details are different: Xu et al., "Seq2seq Fingerprint: An Unsupervised Deep Molecular Embedding for Drug Discovery" (https://doi.org/10.1145/3107411.3107424).
Many types of models require their inputs to have a fixed shape. Since molecules can vary widely in the numbers of atoms and bonds they contain, this makes it hard to apply those models to them. We need a way of generating a fixed length "fingerprint" for each molecule. Various ways of doing this have been designed, such as Extended-Connectivity Fingerprints (ECFPs). But in this example, instead of designing a fingerprint by hand, we will let a `SeqToSeq` model learn its own method of creating fingerprints.
A `SeqToSeq` model performs sequence to sequence translation. For example, they are often used to translate text from one language to another. It consists of two parts called the "encoder" and "decoder". The encoder is a stack of recurrent layers. The input sequence is fed into it, one token at a time, and it generates a fixed length vector called the "embedding vector". The decoder is another stack of recurrent layers that performs the inverse operation: it takes the embedding vector as input, and generates the output sequence. By training it on appropriately chosen input/output pairs, you can create a model that performs many sorts of transformations.
In this case, we will use SMILES strings describing molecules as the input sequences. We will train the model as an autoencoder, so it tries to make the output sequences identical to the input sequences. For that to work, the encoder must create embedding vectors that contain all information from the original sequence. That's exactly what we want in a fingerprint, so perhaps those embedding vectors will then be useful as a way to represent molecules in other models!
Let's start by loading the data. We will use the MUV dataset. It includes 74,501 molecules in the training set, and 9313 molecules in the validation set, so it gives us plenty of SMILES strings to work with.
%% Cell type:code id: tags:
``` python
importdeepchemasdc
tasks,datasets,transformers=dc.molnet.load_muv()
train_dataset,valid_dataset,test_dataset=datasets
train_smiles=train_dataset.ids
valid_smiles=valid_dataset.ids
```
%% Output
About to load MUV dataset.
Loading dataset from disk.
Loading dataset from disk.
Loading dataset from disk.
%% Cell type:markdown id: tags:
We need to define the "alphabet" for our `SeqToSeq` model, the list of all tokens that can appear in sequences. (It's also possible for input and output sequences to have different alphabets, but since we're training it as an autoencoder, they're identical in this case.) Make a list of every character that appears in any training sequence.
%% Cell type:code id: tags:
``` python
tokens=set()
forsintrain_smiles:
tokens=tokens.union(set(cforcins))
tokens=sorted(list(tokens))
```
%% Cell type:markdown id: tags:
Create the model and define the optimization method to use. In this case, learning works much better if we gradually decrease the learning rate. We use an `ExponentialDecay` to multiply the learning rate by 0.9 after each epoch.
Let's train it! The input to `fit_sequences()` is a generator that produces input/output pairs. On a good GPU, this should take a few hours or less.
%% Cell type:code id: tags:
``` python
defgenerate_sequences(epochs):
foriinrange(epochs):
forsintrain_smiles:
yield (s,s)
model.fit_sequences(generate_sequences(40))
```
%% Output
Ending global_step 999: Average loss 72.0029
Ending global_step 1999: Average loss 40.7221
Ending global_step 2999: Average loss 31.5364
Ending global_step 3999: Average loss 26.4576
Ending global_step 4999: Average loss 22.814
Ending global_step 5999: Average loss 19.5248
Ending global_step 6999: Average loss 16.4594
Ending global_step 7999: Average loss 18.8898
Ending global_step 8999: Average loss 13.476
Ending global_step 9999: Average loss 11.5528
Ending global_step 10999: Average loss 10.1594
Ending global_step 11999: Average loss 10.6434
Ending global_step 12999: Average loss 6.57057
Ending global_step 13999: Average loss 6.46177
Ending global_step 14999: Average loss 7.53559
Ending global_step 15999: Average loss 4.95809
Ending global_step 16999: Average loss 4.35039
Ending global_step 17999: Average loss 3.39137
Ending global_step 18999: Average loss 3.5216
Ending global_step 19999: Average loss 3.08579
Ending global_step 20999: Average loss 2.80738
Ending global_step 21999: Average loss 2.92217
Ending global_step 22999: Average loss 2.51032
Ending global_step 23999: Average loss 1.86265
Ending global_step 24999: Average loss 1.67088
Ending global_step 25999: Average loss 1.87016
Ending global_step 26999: Average loss 1.61166
Ending global_step 27999: Average loss 1.40708
Ending global_step 28999: Average loss 1.4488
Ending global_step 29801: Average loss 1.33917
TIMING: model fitting took 5619.924 s
%% Cell type:markdown id: tags:
Let's see how well it works as an autoencoder. We'll run the first 500 molecules from the validation set through it, and see how many of them are exactly reproduced.
Now we'll trying using the encoder as a way to generate molecular fingerprints. We compute the embedding vectors for all molecules in the training and validation datasets, and create new datasets that have those as their feature vectors. The amount of data is small enough that we can just store everything in memory.
Computationally predicting molecular solubility through is useful for drug-discovery. In this tutorial, we will use the `deepchem` library to fit a simple statistical model that predicts the solubility of drug-like compounds. The process of fitting this model involves four steps:
1. Loading a chemical dataset, consisting of a series of compounds along with aqueous solubility measurements.
2. Transforming each compound into a feature vector $v \in \mathbb{R}^n$ comprehensible to statistical learning methods.
3. Fitting a simple model that maps feature vectors to estimates of aqueous solubility.
4. Visualizing the results.
%% Cell type:markdown id: tags:
We need to load a dataset of estimated aqueous solubility measurements [1] into deepchem. The data is in CSV format and contains SMILES strings, predicted aqueaous solubilities, and a number of extraneous (for our purposes) molecular properties. Here is an example line from the dataset:
|Compound ID|ESOL predicted log solubility in mols per litre|Minimum Degree|Molecular Weight|Number of H-Bond Donors|Number of Rings|Number of Rotatable Bonds|Polar Surface Area|measured log solubility in mols per litre|smiles|
Most of these fields are not useful for our purposes. The two fields that we will need are the "smiles" field and the "measured log solubility in mols per litre". The "smiles" field holds a SMILES string [2] that specifies the compound in question. Before we load this data into deepchem, we will load the data into python and do some simple preliminary analysis to gain some intuition for the dataset.
print("Columns of dataset: %s"%str(dataset.columns.values))
print("Number of examples in dataset: %s"%str(dataset.shape[0]))
```
%% Output
ERROR:root:Line magic function `%autoreload` not found.
Automatic pdb calling has been turned OFF
/home/leswing/anaconda3/envs/deepchem27/lib/python2.7/site-packages/sklearn/cross_validation.py:44: DeprecationWarning: This module was deprecated in version 0.18 in favor of the model_selection module into which all the refactored classes and functions are moved. Also note that the interface of the new CV iterators are different from that of this module. This module will be removed in 0.20.
"This module will be removed in 0.20.", DeprecationWarning)
Columns of dataset: ['Compound ID' 'ESOL predicted log solubility in mols per litre'
'Minimum Degree' 'Molecular Weight' 'Number of H-Bond Donors'
'Number of Rings' 'Number of Rotatable Bonds' 'Polar Surface Area'
'measured log solubility in mols per litre' 'smiles']
Number of examples in dataset: 1128
%% Cell type:markdown id: tags:
To gain a visual understanding of compounds in our dataset, let's draw them using rdkit. We define a couple of helper functions to get started.
plt.xlabel('Measured log-solubility in mols/liter')
plt.ylabel('Number of compounds')
plt.title(r'Histogram of solubilities')
plt.grid(True)
plt.show()
```
%% Output
%% Cell type:markdown id: tags:
With our preliminary analysis completed, we return to the original goal of constructing a predictive statistical model of molecular solubility using `deepchem`. The first step in creating such a molecule is translating each compound into a vectorial format that can be understood by statistical learning techniques. This process is commonly called featurization. `deepchem` packages a number of commonly used featurization for user convenience. In this tutorial, we will use ECPF4 fingeprints [3].
`deepchem` offers an object-oriented API for featurization. To get started with featurization, we first construct a ```Featurizer``` object. `deepchem` provides the ```CircularFingeprint``` class (a subclass of ```Featurizer``` that performs ECFP4 featurization).
%% Cell type:code id: tags:
``` python
importdeepchemasdc
featurizer=dc.feat.CircularFingerprint(size=1024)
```
%% Cell type:markdown id: tags:
Now, let's perform the actual featurization. `deepchem` provides the ```CSVLoader``` class for this purpose. The ```featurize()``` method for this class loads data from disk and uses provided ```Featurizer```instances to transform the provided data into feature vectors.
To perform machine learning upon these datasets, we need to convert the samples into datasets suitable for machine-learning (that is, into data matrix $X \in \mathbb{R}^{n\times d}$ where $n$ is the number of samples and $d$ the dimensionality of the feature vector, and into label vector $y \in \mathbb{R}^n$). `deepchem` provides the `Dataset` class to facilitate this transformation. This style lends itself easily to validation-set hyperparameter searches, which we illustate below.
%% Cell type:code id: tags:
``` python
loader=dc.data.CSVLoader(
tasks=["measured log solubility in mols per litre"],smiles_field="smiles",
featurizer=featurizer)
dataset=loader.featurize(dataset_file)
```
%% Output
Loading raw samples now.
shard_size: 8192
About to start loading CSV from ../../datasets/delaney-processed.csv
Loading shard 1 of size 8192.
Featurizing sample 0
Featurizing sample 1000
TIMING: featurizing shard 0 took 2.016 s
TIMING: dataset construction took 2.061 s
Loading dataset from disk.
/home/leswing/Documents/deepchem/deepchem/data/data_loader.py:42: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal
if y[ind, task] == "":
%% Cell type:markdown id: tags:
When constructing statistical models, it's necessary to separate the provided data into train/test subsets. The train subset is used to learn the statistical model, while the test subset is used to evaluate the learned model. In practice, it's often useful to elaborate this split further and perform a train/validation/test split. The validation set is used to perform model selection. Proposed models are evaluated on the validation-set, and the best performed model is at the end tested on the test-set.
Choosing the proper method of performing a train/validation/test split can be challenging. Standard machine learning practice is to perform a random split of the data into train/validation/test, but random splits are not well suited for the purposes of chemical informatics. For our predictive models to be useful, we require them to have predictive power in portions of chemical space beyond the set of molecules in the training data. Consequently, our models should use splits of the data that separate compounds in the training set from those in the validation and test-sets. We use Bemis-Murcko scaffolds [5] to perform this separation (all compounds that share an underlying molecular scaffold will be placed into the same split in the train/test/validation split).
Let's visually inspect some of the molecules in the separate splits to verify that they appear structurally dissimilar. The `FeaturizedSamples` class provides an `itersamples` method that lets us obtain the underlying compounds in each split.
Notice the visual distinction between the train/validation splits. The most-common scaffolds are reserved for the train split, with the rarer scaffolds allotted to validation/test.
%% Cell type:markdown id: tags:
The performance of common machine-learning algorithms can be very sensitive to preprocessing of the data. One common transformation applied to data is to normalize it to have zero-mean and unit-standard-deviation. We will apply this transformation to the log-solubility (as seen above, the log-solubility ranges from -12 to 2).
The next step after processing the data is to start fitting simple learning models to our data. `deepchem` provides a number of machine-learning model classes.
In particular, `deepchem` provides a convenience class, ```SklearnModel``` that wraps any machine-learning model available in scikit-learn [6]. Consequently, we will start by building a simple random-forest regressor that attempts to predict the log-solubility from our computed ECFP4 features. To train the model, we instantiate the ```SklearnModel``` object, then call the ```fit()``` method on the ```train_dataset``` we constructed above. We then save the model to disk.
We next evaluate the model on the validation set to see its predictive power. `deepchem` provides the `Evaluator` class to facilitate this process. To evaluate the constructed `model` object, create a new `Evaluator` instance and call the `compute_model_performance()` method.
The performance of this basic random-forest model isn't very strong. To construct stronger models, let's attempt to optimize the hyperparameters (choices made in the model-specification) to achieve better performance. For random forests, we can tweak `n_estimators` which controls the number of trees in the forest, and `max_features` which controls the number of features to consider when performing a split. We now build a series of `SklearnModel`s with different choices for `n_estimators` and `max_features` and evaluate performance on the validation set.
Model 8/8, Metric r2_score, Validation set 7: 0.157095
best_validation_score so far: 0.275000
computed_metrics: [0.94421662988139621]
Best hyperparameters: (100, 'sqrt')
train_score: 0.944217
validation_score: 0.275000
%% Cell type:markdown id: tags:
The best model achieves significantly higher $R^2$ on the validation set than the first model we constructed. Now, let's perform the same sort of hyperparameter search, but with a simple deep-network instead.
