Classification of molecules by aqueous solubility¶

Here, we will use machine learning to create a simple model to predict whether a molecule is soluble in water or not. This is a classification task, and we will train (and evaluate) different types of models based on a published dataset.

This notebook is inspired by Pat Walters (https://github.com/PatWalters/chem_tutorial/blob/main/tutorial_03_eda.ipynb).

The data we will be using is from the paper “AqSolDB, a curated reference set of aqueous solubility and 2D descriptors for a diverse set of compounds”, M.C. Sorkun, A. Khetan and S. Er, Scientific Data, 2019, 6, 143, DOI: 10.1038/s41597-019-0151-1

The goal of the authors in developing this dataset and making it available is stated in their abstract:

we generated a new database of compounds, AqSolDB, by merging a total of nine different aqueous solubility datasets, curating the merged data, standardizing and validating the compound representation formats, marking with reliability labels, and providing 2D descriptors of compounds as a Supplementary Resource

and their “usage notes” section:

The availability of the calculated 2D descriptors makes it possible to directly use the data for developing machine learning models.

The authors also made available the code that they used to generate the dataset, from other existing datasets, including curation steps: https://doi.org/10.24433/CO.1992938.v1 This is very good open science practice, allowing reproducible research!

If you are interested in this topic, I encourage you to go read the paper and the associated code. Here, we will directly use their data (solubility, list of molecules and associated molecular descriptors) to train a machine learning algorithm for a classification task.

Content of the dataset¶

We have read a local copy of the data set.

We first examine the content of the data set: how many records (molecules) does it have? What information is present for each? What is the type and range of the various descriptors?

Some data fields describe the molecules and their structure:

• Name is the molecule name (not sure if common or IUPAC name)
• InChI and InChIKey are (see Wikipedia)
• SMILES is the (see Wikipedia)

Then there are 17 molecular descriptors, calculated from the molecular structure by RDKit and listed in the paper (Table 2):

Finally, there is data about the solubility. Solubility is the experimental aqueous solubility value (LogS, actually). And because the dataset was obtained by merging together different sources, some solubility values come from different independent measurements. Occurrences is the number of measurements for each molecule, and SD is the standard deviation between measurements. Finally, each molecular is assigned a reliability group (Group).

• G1 - occurs once in the dataset
• G2 - occurs twice in the dataset, values are considered different (SD > 0.5)
• G3 - occurs twice in the dataset, values are considered consistent (SD <= 0.5)
• G4 - occurs 3+ times in the dataset, values are considered different (SD > 0.5)
• G5 - occurs 3+ times in the dataset, values are considered consistent (SD <= 0.5)

Selection / curation of the data¶

We probably want to restrict our analysis to some of the data, keeping only the records that we consider the most reliable. How many molecules are in each group?

We probably want to exclude molecules in the G2 and G4 groups: their solubility has been measured several times, with quite different values, so those values have high uncertainty. But what about G1?

Most of our dataset is from molecules in the G1 group, where the solubility value is derived from a single measurement. For those, we have no idea of the reliability of the data. On the other hand, if we exclude them, our dataset size goes down from 9564 to 1818. What is best for our machine learning model? More data to train on, or less data but with better accuracy?

There is no universal choice, really. In this tutorial, we will start by keeping only the more accurate data, for now.

Exploring the data¶

We can have a look at the distribution of the solubility values in the data set:

Or visually inspect the data to try and see how different descriptors appear to correlate with the solubility:

You can play around with different descriptors… can you guess which is going to be the most important one?

Setting up a classification problem¶

We currently have as data a series of numerical descriptors for each molecule, and a value of solubility. We could therefore create a regression model, which would predict the value of solubility constant (a continuous variable).

But suppose, for our purposes, we are not interested in the exact value of the solubility constant, but only on whether a molecule is soluble or not in water. This is then a classification model, where the outcome we want to predict is a discrete label or category: soluble or insoluble.

We therefore start by creating these categories. We define a molecule as soluble in water if its solubility is higher than 200 μM.

Having modified the whole dataset df, we need to extract again our restricted dataset of high accuracy data df1

And we can check that the data distribution between soluble and insoluble molecules is not too unbalanced:

We see that there are more soluble molecules than insoluble molecules in the dataset, but both sets contain enough molecules that we are able to use the data for prediction.

We are going to train different types of models for this classification task. In order to evaluate their statistical accuracy, we are going to split the entire data set between a training set and a testing set. Scikit-learn has a function for that: train_test_split.

This is a relatively simple way to measure accuracy, and suitable for this simple example. We could also use more advanced procedures, such as cross-validation.

We will now store separately our training features (the descriptors), training targets (the classes to predict: yes/no), testing features and testing targets.

We start by illustrating the process with very simple data, restricted to only two descriptors: molecular weight and number of aromatic rings.

First model: logistic regression¶

As our first model, we choose a very simple example: a logistic regressor, which is a linear model for classification. (And yes, despite its name, it is a classification model and not a regression model. I agree that it is confusing.)

We create a two-step pipeline for this model:

1. standardization of the features
2. logistic regression

In the regressor, we set the penalty parameter to none, so there is no regularization of the model (no constraint on the coefficients of the model). This makes the logistic regression behave as a purely linear model, and will be simpler for discussion later.

Once the model was trained, we measure its accuracy on testing data:

This accuracy is the percentage of correct predictions: it is ~80% in this case, based on test data.

To what do we compare this accuracy? After all, in a game of “heads or tails”, even guessing randomly has a 50% accuracy. Here, we will compare the accuracy to the performance of a so-called dummy classifier, i.e. a model that would always predict the most common label.

Here, this dummy classifier is a model that predicts all molecules to be soluble (because it is the most frequent category). And because 73% of molecules in our test data set are soluble, that “stupid” model has 73% accuracy. So we see that our logistic regression model performance is a little bit better (but not a lot) than the dummy classifier.

Let's see if a different choice of descriptors would make a difference to the model accuracy: we'll take molar weight and this MolLogP descriptor now (more on that below):

The accuracy has definitely gone up, with those two descriptors it's now at about 90%. We can also check and confirm that training accuracy and test accuracy are comparable, so our model does not appear to be overfitting (which is expected from having only two features).

Because we only have two descriptors, we can actually plot the data in this space:

In this 2D feature space, we can also plot the decision function of our trained model. There is no ready-made function for this, so we create our own below. The details of the code are not directly relevant, so you can skip to the resulting graph…

Let’s discuss that graph for a minute. Why does the decision function for the logistic regression model look like this?

As we have seen before, a linear regression model will predict a continuous target, as a linear function of features. But classification tasks can be carried out with linear models, as well. When the target is a binary outcome, we use (like here) a logistic function to model the probability of each outcome, and this is why this model is known as logistic regression. For each point, the classifier will then predict the outcome that it finds most likely.

We plotted above the decision function, red for soluble and blue for insoluble. The edge of that decision function is the curve separating the two areas, where the probability of each is 1/2. In our logistic regression model, it is a straight line.

We noted earlier that we did not impose any regularization by setting the parameter penalty to None. This allows the fitting to allow coefficients as big (or small) as needed. Since the line is oblique, it means that the model uses a combination of both features. We can find the coefficient of each feature in the equation, which will give us an idea of their weight/importance.

And we see that MolLogP is the most important feature in the classification.

Can we make chemical sense of this? MolLogP was calculated by RDKit, where it is described as an “atom-based calculation of the octanol/water partition coefficient log(P)”. More accurately, it is an approximation of log(P), based on this paper by Wildman and Crippen. The approximation is designed to sum contribution from different types of atoms present in a molecule, and coefficients for it were trained against existing experimental data.

Of course, water solubility and octanol/water separation are related, so it makes a lot of sense that MolLogP would be an important descriptor for water solubility.

More descriptors¶

As seen above, the MolLogP molecular descriptor is very well-suited for the prediction of solubility. But using it feels like cheating a bit, since its formula was itself trained by statistical methods on experimental data (which is, in itself, a learning process). So we ask the question: can we achieve the same model accuracy by using other descriptors?

We are making a function that will train a model for each descriptor, and give us as result the accuracy measured on the test dataset:

Many of these descriptors, alone, produce models that are not statistically better than the “dummy classifier” which is our baseline. But how do they perform when put all together?

So we can indeed create a model with descriptors other than MolLogP, that can statistically reach the same accuracy (or almost)!

What are the features that contribute the most to this model?

Let's display this list of coefficients more nicely:

What type of errors are we making? Confusion matrix¶

We have calculated before a statistical accuracy of the predictor. But in a classification task, there are different types of mispredictions, and we may want to check those. For example, in some applications we could be interested to minimize specifically false positives or false negatives.

We can obtain detailed information by calculating and displaying the confusion matrix for the model.

We can see that, in this case, we tend to predict more insoluble molecules as soluble, than we predict soluble molecules as insoluble. We can link this trend to the nature of the training data set, which is biased towards soluble molecules.

Multiclass classification¶

The characterization of confusion matrices is also very interesting when we have multiclass classification tasks. The logistic regression algorithm can also be used for multi-class classification. We are going to test this by defining the solubility of molecules in 3 classes:

• highly soluble (> 1 mM)
• medium solubility (30 μM to 1 mM)
• low solubility (< 30 μM)

We are also going to use now a bigger dataset, this time by introducing the data in G1 group (where only one experimental measurement of solubility is available):

Let’s see how this new data looks like:

The workflow now should be familiar to you, but you can make changes and experiment!

⚠️ That error is not usual, though. What is happening? Because we have more data (and many descriptors), the numerical procedure for the fitting of the linear model did not converge within the default number of iterations. We have to increase that number (max_iter) from 100 to 500:

So the accuracy for our “simple” linear model is now lower for 3-way classification than for binary classification. But what about the confusion matrix?

We can see what is happening: the model has a hard time predicting “medium” solubility: predictions for both low and high solubility have higher accuracy, but most of the records whose true label is medium, are predicted as either high or low (but not medium).

Maybe a linear model is not ideal for this task. Could another type of model perform better?

Using a nearest neighbors model¶

Another category of classification models is that of nearest neighbors methods. It is a type of non-parametric learning, in that it does not attempt to find ideal parameters in an general internal model (like a linear model). Instead, the “training” simply stores a copy of the training data, and classification is done by looking this data up.

Nonparametric methods are good when you have a lot of data and no prior knowledge, and when you don’t want to worry too much about choosing just the right features. (Artificial Intelligence: A Modern Approach)

In neighbors-based methods, classification is computed from a simple majority vote of the nearest neighbors of the point to be predicted. The $k$-nearest neighbors classification (KNeighborsClassifier) is the most common technique: for each new point, it will predict based on a fixed number ($k$) of neighbors in the training data.

We see that, even though the model chosen is very different in nature (kNN instead of logistic regression), the workflow remains the same, and the scikit-learn functions are the same:

• create a pipeline, combining normalization and a ML model
• call the fit function on training data
• call the score function on test data

The accuracy is actually well improved! The $k$-nearest neighbors classifier has an accuracy of 76% (up from 70% for the logistic regression). We see in particular that it improved the quality of prediction for the “medium” solubility significantly… although it is still not ideal.

Going further…¶

Based on the dataset given here, we have introduced classification tasks (both binary and multi-class) and see two algorithms at play. Here are some directions that you might want to explore in the future:

• The $k$-nearest neighbors method depends on one crucial hyperparameter: $k$. How does the number of neighbors impact the model? What happens for very small and very large values of $k$? For this data set, what is the value that provides the best accuracy?
• Another family of classification algorithms are decision trees. They are relatively easy to understand and interpret. Can you train a decision tree for this task? How does the accuracy compare with other models? Does an ensemble method (like a random forest) perform better than a single decision tree?
• How do the models produced for the “small” (highly accurate) data set and “big” data set compare? Can we make statistical comparisons between the values present in both data sets?
• We have set up our problem as a classification task by defining cut-off values of solubility for soluble and insoluble molecules. Could we, as an alternative, use the same solubility data to create a regression task? How would the resulting predictions compare?