QSAR_modelling.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Workshop ULLA2023: QSAR modelling\n",
    "\n",
    "In this workshop you will design your own compounds and predict their properties using QSAR (Quantitative Structure-Activity Relationship) modelling. You will learn how to build and evaluate QSAR models and how to use them to predict the properties of new compounds."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Install QSPRpred"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " <a href=\"https://github.com/CDDLeiden/QSPRPred\">QSPRPred</a> is an open-source software libary for building **Quantitative Structure Property Relationship (QSPR)** models developed by Gerard van Westens Computational Drug Discovery group. It provides a unified interface for building QSPR models based on different types of descriptors and machine learning algorithms.\n",
    " In this tutorial, you will use QSPRPred to build QSAR models for your chosen target protein. You can install QSPRPred using the following command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install git+https://github.com/CDDLeiden/QSPRPred.git@BOO-2023"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create your dataset"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, you need to create a dataset for your target protein. The data will be collected from <a href=\"https://github.com/CDDLeiden/Papyrus-scripts\">Papyrus</a>. \n",
    "\n",
    "The Papyrus dataset is a highly curated compilation of bioactivity data intended for modeling in drug discovery. Apart from the bioactivity data contained in the [ChEMBL database](https://www.ebi.ac.uk/chembl/), the Papyrus dataset contains binary data for classification tasks from the [ExCAPE-DB](https://solr.ideaconsult.net/search/excape/), and bioactivity data from a number of kinase-specific papers. The Papyrus dataset consists of almost 60M compound-protein pairs, representing data of around 1.2M unique compounds and 7K proteins across 499 different organisms.\n",
    "\n",
    "The aggregated bioactivity data is standardized, repaired, and normalized to form the Papyrus dataset. The Papyrus dataset contains \"high quality\" data associated with pChEMBL values that represent exact bioactivity values measured and associated with a single protein or complex subunit. A pChEMBL value is a canonical activity metric defined as $-log_{10}(molar IC_{50}, XC_{50}, EC_{50}, AC_{50}, Ki, Kd, or potency)$. Moreover, \"low and medium quality\" data that is associated with potentially multiple proteins or a homologous single proteins and/or only approximate bioacivity values or binary active/inactive label is also included in the Papyrus dataset, but is not used in this workshop.\n",
    "\n",
    "Here I've already made a pre-selection from the Papyrus dataset (including only the data for all the proteins used in the workshop) to shorten the data collection process. You still need to collect the data from the dataset for your specific target protein. You can do this by replacing the accession number in the code below with your target protein accession number. You can find the accession number of your target protein in the <a href=\"https://www.uniprot.org/\">UniProt</a> database.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in all data\n",
    "\n",
    "import pandas as pd\n",
    "df = pd.read_csv('data/papyrus_data.tsv', sep='\\t')\n",
    "\n",
    "# filter data for target of interest\n",
    "MY_TARGET = 'P21802' # REPLACE WITH YOUR TARGET ACCESSION\n",
    "\n",
    "df = df[df['accession'] == MY_TARGET]\n",
    "\n",
    "# keep only high quality data\n",
    "df = df[df['Quality'].isin(['High'])] # if you want to keep also medium and low quality data, remove this line\n",
    "\n",
    "# Create molecule table for visualization\n",
    "from qsprpred.data.data import MoleculeTable\n",
    "\n",
    "mt = MoleculeTable(df=df, name=MY_TARGET, store_dir='data')\n",
    "\n",
    "mt.getDF()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparing data for modelling"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After you have collected the data for your target protein, you need to prepare the data for modelling.\n",
    "First we will specify the name of the target property (pchembl_value_Median) and the model task (REGRESSION) to \n",
    "create our QSPRdataset object. The QSPRdataset object is a wrapper around the pandas dataframe that contains the data for modelling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qsprpred.models.tasks import TargetTasks\n",
    "from qsprpred.data.data import QSPRDataset\n",
    "\n",
    "target_props=[{\n",
    "                \"name\": \"pchembl_value_Median\", # name of the target column in the dataset\n",
    "                \"task\": TargetTasks.REGRESSION, # specify the task type (SINGLECLASS, MULTICLASS, REGRESSION)\n",
    "                }]\n",
    "\n",
    "# Create a QSPRDataset instance used for training and evaluation of QSPR models\n",
    "dataset = QSPRDataset.fromMolTable(mt, target_props=target_props)\n",
    "dataset.targetProperties"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we will use the `prepareDataset` method to prepare the data for modelling. This method will perform the following steps:\n",
    "\n",
    "1. Standardize the SMILES sequences\n",
    "2. Split the data into training and test sets (80/20 split)\n",
    "3. Calculate the descriptors (Morgan Fingerprints) as discussed in the lecture (See also image below)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![descriptors](figures/descriptors.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qsprpred.data.utils.descriptorsets import FingerprintSet\n",
    "from qsprpred.data.utils.descriptorcalculator import MoleculeDescriptorsCalculator\n",
    "from qsprpred.data.utils.datasplitters import randomsplit\n",
    "\n",
    "# Calculate MorganFP\n",
    "feature_calculator = MoleculeDescriptorsCalculator(descsets = [FingerprintSet(fingerprint_type=\"MorganFP\", radius=3, nBits=2048)])\n",
    "\n",
    "# Do a random split for creating the train (80%) and test set (20%)\n",
    "rand_split = randomsplit(0.2)\n",
    "\n",
    "# calculate compound features and split dataset into train and test\n",
    "dataset.prepareDataset(\n",
    "    split=rand_split,\n",
    "    feature_calculators=[feature_calculator]\n",
    ")\n",
    "\n",
    "# Let's have a look at the features of the first five compounds\n",
    "print(f\"Number of features (bits in the MorganFP): {len(dataset.getFeatures()[0].columns)}\")\n",
    "display(dataset.getFeatures()[0].head())\n",
    "\n",
    "print(f\"Number of samples train set: {len(dataset.y)}\")\n",
    "print(f\"Number of samples test set: {len(dataset.y_ind)}\")\n",
    "\n",
    "# Let's save the dataset for later\n",
    "dataset.save()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data Visualization"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have prepared the data for modelling, we can visualize the data to get a better understanding of what is in our dataset. We will first plot a histogram of the target property values to get an idea of the distribution of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create histogram of pchembl values in the dataset\n",
    "import seaborn as sns\n",
    "sns.histplot(dataset.getDF()['pchembl_value_Median'], bins=20)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the goal of this workshop is to design your own compounds, we will also have a look at the compounds with the highest affinity for your target. We will use the `draw.MolsToGridImage` method to draw the compounds with the highest pChEMBL values. You can change the number of compounds to draw by changing the `NUM_COMPOUNDS` variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize the compounds with the highest pchembl values in the dataset\n",
    "from rdkit.Chem import Draw\n",
    "from rdkit.Chem.Draw import IPythonConsole\n",
    "from rdkit.Chem import MolFromSmiles\n",
    "\n",
    "# number of compounds to show\n",
    "NUM_COMPOUNDS = 10\n",
    "\n",
    "# Sort the dataset by pchembl value\n",
    "dataset_sorted = dataset.getDF().sort_values(by='pchembl_value_Median', ascending=False)\n",
    " \n",
    "# show average pchembl value per scaffold and the count of compounds per scaffold\n",
    "Draw.MolsToGridImage([MolFromSmiles(smiles) for smiles in dataset_sorted[:NUM_COMPOUNDS].SMILES], molsPerRow=5, subImgSize=(200,200), legends=[f\"{dataset_sorted['pchembl_value_Median'][idx]:.2f}\" for idx in dataset_sorted[:NUM_COMPOUNDS].index])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although molecules with a high pChEMBL value can be a good starting point for designing your own compound, it can be difficult to know what to change to improve the affinity. In drug discovery, chemists often explore the structure-activity relationship (SAR) of a compound by making small changes to a core structure (`scaffold`), such as replacing a substituent or adding a functional group. We can use this information by extracting the scaffolds from the compounds and finding the scaffolds with on average the highests pChEMBL values to better understand the SAR of the compounds in our dataset. \n",
    "\n",
    "Here we will visualize those scaffolds. You can change the number of scaffolds to draw by changing the `NUM_SCAFFOLDS` variable and the `MIN_COMPOUNDS` variable to change the minimum number of compounds that should contain the scaffold. Note, the numbers under the scaffolds are \n",
    "*{index}: {average pChEMBL value} ({number of compounds that contain the scaffold})*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qsprpred.data.utils.scaffolds import Murcko, BemisMurcko\n",
    "\n",
    "# Show top n scaffolds with at least x compounds\n",
    "NUM_SCAFFOLDS = 10\n",
    "MIN_COMPOUNDS = 5\n",
    "\n",
    "dataset.addScaffolds([Murcko()])\n",
    "\n",
    "# get average pchembl value per scaffold\n",
    "scaffolds = dataset.getDF().groupby('Scaffold_Murcko')['pchembl_value_Median'].mean().sort_values(ascending=False)\n",
    "scaffolds = scaffolds.rename('Average pchembl value')\n",
    "\n",
    "# add the number of compounds per scaffold\n",
    "scaffolds = pd.concat([scaffolds, dataset.getDF().groupby('Scaffold_Murcko')['pchembl_value_Median'].count()], axis=1)\n",
    "scaffolds = scaffolds.rename(columns={'pchembl_value_Median': 'Count'})\n",
    "\n",
    "# Drop scaffolds with less than MIN_COMPOUNDS compounds\n",
    "scaffolds = scaffolds[scaffolds['Count'] > MIN_COMPOUNDS]\n",
    " \n",
    "# show average pchembl value per scaffold and the count of compounds per scaffold\n",
    "Draw.MolsToGridImage([MolFromSmiles(scaffold) for scaffold in scaffolds.index[:NUM_SCAFFOLDS]], molsPerRow=5, subImgSize=(200,200), legends=[f\"{idx}: {scaffolds['Average pchembl value'][scaffold]:.2f} ({scaffolds['Count'][scaffold]})\" for idx, scaffold in enumerate(scaffolds.index[:NUM_SCAFFOLDS])])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After you have studied the promising scaffolds, you can select one of the scaffolds to design your own compound. Let's have a look at what molecules that have been tested that contain this scaffold. You can do this by changing the `SCAFFOLD_INDEX` variable to the scaffold that you want to see the compounds of. The index of the scaffold is the first number under the scaffold image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show all compounds from from scaffold N.\n",
    "\n",
    "SCAFFOLD_INDEX = 0 # REPLACE WITH YOUR SCAFFOLD INDEX\n",
    "\n",
    "scaffold = scaffolds.index[SCAFFOLD_INDEX]\n",
    "\n",
    "# get all compounds from scaffold\n",
    "scaffold_df = dataset.getDF()[dataset.getDF()['Scaffold_Murcko'] == scaffold]\n",
    "\n",
    "# sort compounds by pchembl value\n",
    "scaffold_df = scaffold_df.sort_values(by='pchembl_value_Median', ascending=False)\n",
    "\n",
    "# visualize compounds\n",
    "Draw.MolsToGridImage([MolFromSmiles(smiles) for smiles in scaffold_df.SMILES], molsPerRow=5, subImgSize=(200,200), legends=[f\"{scaffold_df['pchembl_value_Median'][idx]:.2f}\" for idx in scaffold_df.index])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Training a ML model"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After preparing the data for modelling, we can train a machine learning model to predict the target property.\n",
    "\n",
    "\n",
    "First, we will specify the machine learning algorithm that we want to use to train our model. For this workshop we will use the Partial Least Squares Regression model from <a href=\"https://scikit-learn.org/stable/\">scikit-learn</a> as this model is fast to train. However, the model algorithm could easily be changed to a different model by replacing the alg argument with a different model from scikit-learn. You can find a list of models in the <a href=\"https://scikit-learn.org/stable/supervised_learning.html\">scikit-learn documentation</a>.\n",
    "\n",
    "See the image below for an overview of the steps involved in training a machine learning model.\n",
    "1. We will first optimize the hyperparameters of the model using the bayesian optimization (See <a href=\"https://towardsdatascience.com/a-conceptual-explanation-of-bayesian-model-based-hyperparameter-optimization-for-machine-learning-b8172278050f\"> this tutorial </a> for an explanation) for five trials. Normally, you would use more trials to find the best hyperparameters, but this can take a long time, so we will limit the number of trials for this workshop. You can change the number of trials by changing the `n_trials` variable.\n",
    "2. Using the `evaluate` function, we will train the model using the optimized hyperparameters on the training data using cross-validation.\n",
    "3. Also in the `evaluate` function, we will evaluate the model on the test data using the optimized hyperparameters.\n",
    "4. Finally we will train the model on all the data with the `fit` method to use it for predicting the activities of new compounds.\n",
    "\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![modelling](figures/modelling_workflow.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qsprpred.models.models import QSPRsklearn\n",
    "from sklearn.cross_decomposition import PLSRegression\n",
    "from qsprpred.models.hyperparam_optimization import OptunaOptimization\n",
    "\n",
    "# This is an SKlearn model, so we will initialize it with the QSPRsklearn class\n",
    "model = QSPRsklearn(base_dir = '.', data=dataset, alg = PLSRegression, name='PLS_REG')\n",
    "\n",
    "# We will first optimize the hyperparameters (n_components and scale) through bayes optimization\n",
    "# the best hyperparameter combination will be saved in PLS_REG_params.json\n",
    "search_space_bs = {\"n_components\": [\"int\", 1, 30], \"scale\": [\"categorical\", [True, False]]}\n",
    "bayesoptimizer = OptunaOptimization(scoring = model.score_func, param_grid=search_space_bs, n_trials=5)\n",
    "best_params = bayesoptimizer.optimize(model)\n",
    "\n",
    "#Then we will evaluate the performance of the best model using the independent test set\n",
    "_ = model.evaluate()\n",
    "\n",
    "# Finally, we need to fit the model on the complete dataset if we want to use it further\n",
    "# model is saved under qspr/models/PLS_REG.json\n",
    "model.fit()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After training a model, you can visualize the predictions on the evaluation folds of the cross validation and the test set to get an impression of the performance of the model by creating a scatter plot of the predicted vs. the actual values. Also the R2 and RMSE values are printed to get a better idea of the performance of the model (cv stands for cross-validation and ind for the test set).\n",
    "\n",
    "- **Coefficient of determination (**$R^{2}$ **score)**: Represents the portion of variance of the target variable that has been explained by the independent variables (features) in the model. $R^{2}$ score varies between 1.0 (best score) and minus infinite, where 0.0 represents a model that always predicts the average target variable. As the variance is dataset-dependent, it might not be a meaningful metric to compare between datasets. When dealing with linear regression, and model fitting and evaluation are performed on a single dataset, $R^{2}$ is equivalent to the square of the Pearson correlation coefficient and can be noted as $r^{2}$.,\n",
    "- **Root mean square error (RMSE)**: It is the square root of the MSE and represents the standard deviation of the prediction errors with respect to the line of best fit. RMSE is a measure of accuracy and it cannot be applied to compare  between datasets, as it is scale-dependent. It varies between infinite and 0.0 (best).,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qsprpred.plotting.regression import CorrelationPlot\n",
    "\n",
    "plt = CorrelationPlot([model])\n",
    "axes, summary = plt.make(save=False, property_name='pchembl_value_Median')\n",
    "axes[0]\n",
    "\n",
    "print(summary)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Make predictions for your own compounds"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that you have trained a model, you can use it to predict the activities of your own compounds. You can do this by replacing the SMILES sequences in the `smiles` list with your own SMILES sequences. To create your own compounds you can draw them using <a href=https://molview.org>MolView</a> and go to the `Tools` tab and click on `Information card`. You can then copy the SMILES sequence from the `Information card` and paste it in the `list_of_smiles`. The model prediction will then be printed for each compound."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# replace with your own compounds\n",
    "list_of_smiles = ['OCCc1ccn2cnccc12',\n",
    "                  'C1CC1Oc1cc2ccncn2c1',\n",
    "                  'CNC(=O)c1nccc2cccn12'] # REPLACE WITH YOUR OWN COMPOUNDS\n",
    "\n",
    "# make predictions with the model\n",
    "predictions = model.predictMols(list_of_smiles)\n",
    "\n",
    "# show molecules with predicted values using rdkit\n",
    "from rdkit import Chem\n",
    "from rdkit.Chem import Draw\n",
    " \n",
    "mols = [Chem.MolFromSmiles(smi) for smi in list_of_smiles]\n",
    "Draw.MolsToGridImage(mols, molsPerRow=4, subImgSize=(200, 200), legends=[f'{pred[0]:.3f}' for pred in predictions])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ulla2023",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.3"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}