02_molecular_property_prediction_solution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "857f895e4c2ded45",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "# Molecular property prediction"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47292e9e-f8ec-41e3-b685-31e861a35d9f",
   "metadata": {},
   "source": [
    "Install dependencies if necessary (uncomment and run):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2de91fa2dc3b3443",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# !pip install -r requirements.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9923fb4b525f8889",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## 1. Molecular property prediction dataset\n",
    "\n",
    "We will use BACE dataset from the [MoleculeNet benchmark](https://arxiv.org/abs/1703.00564). Throughout this notebook, we will also make heavy use of [scikit-fingerprints library](https://github.com/scikit-fingerprints/scikit-fingerprints), which is a scikit-learn compatible library built around RDKit.\n",
    "\n",
    "The task is classifying inhibitors of Beta-Secretase 1 - a protein enzyme playing a significant role in development of Alzheimer’s disease. This is binary graph classification: molecule inhibits protein production or not.\n",
    "\n",
    "For more information, see: [\"Computational Modeling of β-Secretase 1 (BACE-1) Inhibitors Using Ligand Based Approaches\" G. Subramanian et al.](https://pubs.acs.org/doi/10.1021/acs.jcim.6b00290)\n",
    "\n",
    "scikit-fingerprints has its own data loaders for several popular datasets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7a730860462b9a3",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from skfp.datasets.moleculenet import load_bace\n",
    "\n",
    "smiles_list, y = load_bace()\n",
    "\n",
    "print(f\"Example molecule: {smiles_list[0]}\")\n",
    "print(f\"Example class: {y[0]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad27702a3ef5a42b",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## 2. Preprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c867a4f-1cf2-46f0-8372-0d34499d01d1",
   "metadata": {},
   "source": [
    "### SMILES -> Mol conversion\n",
    "\n",
    "We already covered converting SMILES to RDKit `Mol` objects with pure RDKit. Its function works on one string at a time.\n",
    "\n",
    "scikit-fingerprints allows us to do that for entire dataset with `MolFromSmilesTransformer`. It has the same `.transform()` method as other scikit-learn objects. Note that most classes in scikit-fingerprints do not need `.fit()` call before first usage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a026e271920bafce",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from skfp.preprocessing import MolFromSmilesTransformer\n",
    "\n",
    "mol_from_smiles = MolFromSmilesTransformer()\n",
    "\n",
    "mols = mol_from_smiles.transform(smiles_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eeedb1d3-0cd7-4d80-862d-09b798a7c76f",
   "metadata": {},
   "source": [
    "### Scaffold Split\n",
    "\n",
    "In molecular property prediction, we typically **don't** use random or stratified random split.\n",
    "\n",
    "In the real-world drug design problems, a trained ML model has to perform well on newly designed molecules. They differ significantly from the ones seen in the training set, e.g. to be patentable. We need a splitting strategy that will force and test **out-of-distribution (OOD) generalization**.\n",
    "\n",
    "We can group the molecules by the similarity of their core internal structure, called **scaffold**. This effectively splits the data into sets that differ from one another.\n",
    "\n",
    "It's important to note that there are some slight variations of scaffold split. They differ in definition what is a \"core\" part of the molecule. Most benchmarks provide explicit splits for datasets, e.g. Open Graph Benchmark ([OGB](https://ogb.stanford.edu/)) provides standardized scaffold splits for MoleculeNet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8beacf48077f097c",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from skfp.model_selection import scaffold_train_test_split\n",
    "\n",
    "\n",
    "train_idxs, test_idxs = scaffold_train_test_split(\n",
    "    mols, test_size=0.2, return_indices=True\n",
    ")\n",
    "\n",
    "# split mols and labels\n",
    "mols_train = np.array(mols)[train_idxs]\n",
    "mols_test = np.array(mols)[test_idxs]\n",
    "\n",
    "y_train = y[train_idxs]\n",
    "y_test = y[test_idxs]\n",
    "\n",
    "print(f\"Train set size: {len(mols_train)}\")\n",
    "print(f\"Test set size: {len(mols_test)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "511d9323ab7e81a0",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "mols_train[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a05d04588835999a",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from rdkit.Chem.Scaffolds.MurckoScaffold import GetScaffoldForMol\n",
    "\n",
    "\n",
    "GetScaffoldForMol(mols_train[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2a7fba3-7867-4019-928e-052c27e12886",
   "metadata": {},
   "source": [
    "### Standardize\n",
    "\n",
    "Molecular standardizer performs basic sanitization and standardization of the molecule object. It will make sure that e.g. certain functional groups or electric charges are represented in a uniform way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "677c1711bc2942fe",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from skfp.preprocessing import MolStandardizer\n",
    "\n",
    "\n",
    "standardizer = MolStandardizer()\n",
    "\n",
    "mols_train = standardizer.transform(mols_train)\n",
    "mols_test = standardizer.transform(mols_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8f02792d5f4d050",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## 3. Classification with molecular fingerprints\n",
    "\n",
    "To transform molecules into feature vectors, we will moleculecular fingerprints. scikit-fingeprints implements a lot of those.\n",
    "\n",
    "We will start with popular ECFP fingerprint, which is a hashed fingerprint, using circular substructures. This turns our dataset into typical tabular classification problem. Then, we can use any off-the-shelf classifier, like Random Forest.\n",
    "\n",
    "This is a single-task dataset. In chemistry, we often have **multitask** datasets, where we do many classifications at once. For example, molecule can be toxic in many different ways. `multioutput_auroc_score` in scikit-fingerprints will work in both single-task and multitask cases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "386ab00e-428b-428a-86c3-004ec2e31b84",
   "metadata": {},
   "outputs": [],
   "source": [
    "from skfp.fingerprints import ECFPFingerprint\n",
    "\n",
    "\n",
    "# create fingerprint transformer object\n",
    "ecfp_fp = ECFPFingerprint()\n",
    "\n",
    "# transform molecules into feature vectors\n",
    "X_train_ecfp = ecfp_fp.transform(mols_train)\n",
    "X_test_ecfp = ecfp_fp.transform(mols_test)\n",
    "\n",
    "print(f\"Fingerprint data shape: {X_train_ecfp.shape}\")\n",
    "print(f\"Example vector: {X_train_ecfp[0]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3a33efe-a209-4a1c-a098-4f05becb3a27",
   "metadata": {},
   "source": [
    "Let's train the classifier, using fingerprint features as inputs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48ee7255-238b-4002-a4c0-d7ffb28e61e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from skfp.metrics import multioutput_auroc_score\n",
    "\n",
    "\n",
    "clf = RandomForestClassifier(n_jobs=-1, random_state=0)\n",
    "clf.fit(X_train_ecfp, y_train)\n",
    "\n",
    "y_pred = clf.predict_proba(X_test_ecfp)[:, 1]\n",
    "auroc = multioutput_auroc_score(y_test, y_pred)\n",
    "\n",
    "print(f\"ECFP AUROC for Random Forest: {auroc:.2%}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ab5e102-e35b-4fc0-9452-b3a7d1101d8f",
   "metadata": {},
   "source": [
    "We can also other models, such as k-nearest-neighbours. Note that we have either binary or count vectors, so we should use an appropriate distance. It is typically Tanimoto distance\n",
    "\n",
    "Binary Tanimoto distance formula:\n",
    "\n",
    "$$\n",
    "\\text{dist}(\\vec{a}, \\vec{b}) = 1 - \\frac{|\\vec{a} \\cap \\vec{b}|}{|\\vec{a}| + |\\vec{b}| - |\\vec{a} \\cap \\vec{b}|}\n",
    "$$\n",
    "\n",
    "We can also define Tanimoto distance for count data:\n",
    "\n",
    "$$\n",
    "\\text{dist}(\\vec{a}, \\vec{b}) = 1 - \\frac{\\vec{a} \\cdot \\vec{b}}{\\|\\vec{a}\\|^2 + \\|\\vec{b}\\|^2 - \\vec{a} \\cdot \\vec{b}}\n",
    "$$\n",
    "\n",
    "If you're interested why Tanimoto distance is used, see e.g. [\"Why is Tanimoto index an appropriate choice for fingerprint-based similarity calculations?\" D. Bajusz et al.](https://doi.org/10.1186/s13321-015-0069-3)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d70e7658ecf7a1d1",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "### Exercise 1\n",
    "\n",
    "Implement kNN classifier, using count ECFP fingerprints and Tanimoto distance. Use [scikit-fingerprints documentation on Tanimoto distance](https://scikit-fingerprints.github.io/scikit-fingerprints/) if necessary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c32519b4-1c04-4196-b2d8-de95b0c2f22b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import Tanimoto count distance from scikit-fingerprints\n",
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "from skfp.distances import tanimoto_count_distance\n",
    "\n",
    "# create kNN model with appropriate metric\n",
    "clf = KNeighborsClassifier(n_jobs=-1, metric=tanimoto_count_distance)\n",
    "\n",
    "# fit, predict\n",
    "clf.fit(X_train_ecfp, y_train)\n",
    "y_pred = clf.predict_proba(X_test_ecfp)[:, 1]\n",
    "\n",
    "# calculate and print AUROC score\n",
    "print(f\"ECFP AUROC for kNN: {multioutput_auroc_score(y_test, y_pred):.2%}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "433a208485dd03c9",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "### Exercise 2\n",
    "\n",
    "scikit-fingerprints implements over 30 different algorithms. Another popular choice besides ECFP is MACCS, a substructure-based fingerprint.\n",
    "\n",
    "Implement Random Forest model on binary MACCS fingerprint. Use [scikit-fingerprints documentation](https://scikit-fingerprints.github.io/scikit-fingerprints/) as needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5899f92f47c88477",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# import MACCS fingerprint\n",
    "from skfp.fingerprints import MACCSFingerprint\n",
    "\n",
    "# create MACCS transformer, calculate fingerprints\n",
    "maccs_fp = MACCSFingerprint(n_jobs=-1)\n",
    "\n",
    "X_train_maccs = maccs_fp.transform(mols_train)\n",
    "X_test_maccs = maccs_fp.transform(mols_test)\n",
    "\n",
    "# train Random Forest, predict\n",
    "clf = RandomForestClassifier(n_jobs=-1, random_state=0)\n",
    "clf.fit(X_train_maccs, y_train)\n",
    "\n",
    "y_pred = clf.predict_proba(X_test_maccs)[:, 1]\n",
    "\n",
    "# calculate and print AUROC score\n",
    "print(f\"MACCS AUROC for RF: {multioutput_auroc_score(y_test, y_pred):.2%}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6d20bbc268a2349",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## 4. Hyperparameter tuning"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cb0e882-f0e3-48b6-b1b9-fe3fc77b940e",
   "metadata": {},
   "source": [
    "Molecular fingerprints have many hyperparameters that can be tuned. It's rarely done in papers, but can sometimes greatly improve performance.\n",
    "\n",
    "Examples are:\n",
    "- binary vs count variant\n",
    "- length / number of bits\n",
    "- ECFP radius\n",
    "\n",
    "scikit-fingerprints implements `FingerprintEstimatorGridSearch` class for this. It can speed up the tuning of fingerprint and classifier, by avoiding unnecessary recalculations of the fingerprint.\n",
    "\n",
    "Firstly, we will prepare tuning of the classifier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "233d71c3bdf3e508",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.metrics import make_scorer\n",
    "\n",
    "clf = RandomForestClassifier(n_jobs=-1, random_state=0)\n",
    "\n",
    "clf_param_grid = {\"n_estimators\": [200, 500, 1000]}\n",
    "\n",
    "scorer = make_scorer(multioutput_auroc_score, greater_is_better=True)\n",
    "\n",
    "gridsearch_cv = GridSearchCV(\n",
    "    estimator=clf,\n",
    "    param_grid=clf_param_grid,\n",
    "    scoring=scorer,\n",
    "    verbose=2,\n",
    "    cv=5,\n",
    "    n_jobs=-1,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6331d96c-0bef-47d9-8813-b51cc32a9d88",
   "metadata": {},
   "source": [
    "Now, we will set up the tuning of the fingerprint. `FingerprintEstimatorGridSearch` takes the fingeprints, its grid search, and the classifier grid search tuning objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cde79eeda1e8fb81",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "from skfp.model_selection import FingerprintEstimatorGridSearch\n",
    "\n",
    "ecfp_fp = ECFPFingerprint()\n",
    "\n",
    "fp_grid = {\n",
    "    \"fp_size\": [1024, 2048, 4096],\n",
    "    \"radius\": [2, 3],\n",
    "    \"count\": [False, True],\n",
    "}\n",
    "\n",
    "fp_estimator_cv = FingerprintEstimatorGridSearch(\n",
    "    fingerprint=ecfp_fp,\n",
    "    fp_param_grid=fp_grid,\n",
    "    estimator_cv=gridsearch_cv,\n",
    "    greater_is_better=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe02da2b-789a-4651-9c67-65a9a7ac18ad",
   "metadata": {},
   "source": [
    "Let's perform the tuning now. We will also check the time, best hyperparameters, and performance after tuning."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48a7ffef6e7908ed",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from time import time\n",
    "\n",
    "start = time()\n",
    "fp_estimator_cv.fit(mols_train, y_train)\n",
    "end = time()\n",
    "print(f\"scikit-learn tuning time : {end - start:.2f}s\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac19630a-29b8-41d6-a92a-4e021191a783",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Best fingerprint hyperparameters:\", fp_estimator_cv.best_fp_params_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7afd708-dc0e-4a08-ac70-5480e0d89f48",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Best model hyperparameters:\", fp_estimator_cv.best_estimator_cv_.best_params_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d3c3a1fe4670abf",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "y_pred = fp_estimator_cv.predict_proba(mols_test)[:, 1]\n",
    "\n",
    "print(f\"ECFP AUROC : {multioutput_auroc_score(y_test, y_pred):.2%}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34ccd5e0387844c0",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## 5. Multioutput prediction"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e09d2b3-1080-443c-93d8-55ccdd901a05",
   "metadata": {},
   "source": [
    "Some molecular datasets focus on multiple tasks at a time. This often results in better classifiers, due to more general knowledge and built-in regularization. Random Forest from scikit-learn, combined with molecular fingerprints, is a natural solution for such problems.\n",
    "\n",
    "One of such multi-output datasets is SIDER dataset. Tasks are related to adverse drug reactions (ADRs), or drug side effects, to 27 system organ classes of MedDRA classification.\n",
    "\n",
    "For details, see [\"Low Data Drug Discovery with One-Shot Learning\" H. Altae-Tran et al.](https://pubs.acs.org/doi/10.1021/acscentsci.6b00367)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a479eff06101b1b",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "### Exercise 3\n",
    "\n",
    "Perform model training with molecular fingerprints on SIDER dataset.\n",
    "\n",
    "This consists of steps:\n",
    "1. Load data\n",
    "2. Parse as RDKit molecules\n",
    "3. Scaffold split\n",
    "4. Computing fingerprints\n",
    "5. Training Random Forest\n",
    "6. Evaluation\n",
    "\n",
    "Code templates have been prepared for you below. Note that they sometimes assume variable names. `extract_pos_proba` function will be useful. Use the previous notebook code and [scikit-fingerprints docs](https://scikit-fingerprints.github.io/scikit-fingerprints/) as necessary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d56b694ce00a126",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# import necessary things\n",
    "from skfp.datasets.moleculenet import load_sider\n",
    "\n",
    "# load SIDER dataset: SMILES strings and labels\n",
    "smiles_list_sider, y = load_sider()\n",
    "\n",
    "print(f\"Example molecule: {smiles_list[0]}\")\n",
    "print(f\"Example classes: {y[0]}\")\n",
    "print(f\"Number of outputs: {y[0].shape[0]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf975f9c0208cc85",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# import necessary things\n",
    "from skfp.model_selection import scaffold_train_test_split\n",
    "from skfp.preprocessing import MolFromSmilesTransformer, MolStandardizer\n",
    "\n",
    "# parse SMILES as molecules\n",
    "mol_from_smiles = MolFromSmilesTransformer()\n",
    "mols = mol_from_smiles.transform(smiles_list_sider)\n",
    "\n",
    "# scaffold split with 80-20% proportion\n",
    "train_idxs, test_idxs = scaffold_train_test_split(\n",
    "    mols, test_size=0.2, return_indices=True\n",
    ")\n",
    "\n",
    "mols_train = np.array(mols)[train_idxs]\n",
    "mols_test = np.array(mols)[test_idxs]\n",
    "\n",
    "y_train = y[train_idxs]\n",
    "y_test = y[test_idxs]\n",
    "\n",
    "# create standardizer and standardize molecules\n",
    "standardizer = MolStandardizer()\n",
    "\n",
    "mols_train = standardizer.transform(mols_train)\n",
    "mols_test = standardizer.transform(mols_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "414e53616c0524ed",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# import necessary things\n",
    "from skfp.fingerprints import ECFPFingerprint\n",
    "from skfp.metrics import extract_pos_proba, multioutput_auroc_score\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "\n",
    "# compute fingerprints\n",
    "ecfp_fp = ECFPFingerprint()\n",
    "\n",
    "X_train_ecfp = ecfp_fp.transform(mols_train)\n",
    "X_test_ecfp = ecfp_fp.transform(mols_test)\n",
    "\n",
    "# train RF classifier\n",
    "rf_clf = RandomForestClassifier(n_jobs=-1, random_state=0)\n",
    "rf_clf.fit(X_train_ecfp, y_train)\n",
    "\n",
    "# predict, extract positive class probabilities\n",
    "y_pred = rf_clf.predict_proba(X_test_ecfp)\n",
    "y_pred = extract_pos_proba(y_pred)\n",
    "\n",
    "# calculate AUROC\n",
    "auroc = multioutput_auroc_score(y_test, y_pred)\n",
    "print(f\"SIDER ECFP+RF AUROC: {auroc:.2%}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e80df4d56347bd01",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## 6. Pipelines"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7ac80e1-6a5f-42c4-a223-91e13fc78940",
   "metadata": {},
   "source": [
    "In ML, we often have multi-step pipelines, with preprocessing, training multiple models, merging their ensembles etc. We can do the same with scikit-fingerprints. Using many different fingerprints often helps.\n",
    "\n",
    "Here, we will use 80-10-10% train-valid-test scaffold split, to be able to compare to other papers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9396ab10-6517-47bc-96b1-f645dc8196fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "from skfp.datasets.moleculenet import load_bace\n",
    "from skfp.model_selection import scaffold_train_valid_test_split\n",
    "\n",
    "smiles_list, y = load_bace()\n",
    "\n",
    "train_idxs, valid_idxs, test_idxs = scaffold_train_valid_test_split(\n",
    "    smiles_list, train_size=0.8, valid_size=0.1, test_size=0.1, return_indices=True\n",
    ")\n",
    "\n",
    "# split mols and labels\n",
    "smiles_train = np.array(smiles_list)[train_idxs]\n",
    "smiles_test = np.array(smiles_list)[test_idxs]\n",
    "\n",
    "y_train = y[train_idxs]\n",
    "y_test = y[test_idxs]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f44c27ebd63b6201",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from skfp.fingerprints import ECFPFingerprint, TopologicalTorsionFingerprint\n",
    "from sklearn.pipeline import Pipeline, FeatureUnion\n",
    "from sklearn.preprocessing import MinMaxScaler\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "\n",
    "fps_union = FeatureUnion(\n",
    "    [\n",
    "        (\"tt_fp\", TopologicalTorsionFingerprint()),\n",
    "        (\"ecfp_fp\", ECFPFingerprint()),\n",
    "    ]\n",
    ")\n",
    "\n",
    "rf_clf = RandomForestClassifier(class_weight=\"balanced\", n_jobs=-1, random_state=0)\n",
    "\n",
    "pipeline = Pipeline(\n",
    "    [\n",
    "        (\"mol_from_smiles\", MolFromSmilesTransformer()),\n",
    "        (\"mol_standardizer\", MolStandardizer()),\n",
    "        (\"fps_union\", fps_union),\n",
    "        (\"scaler\", MinMaxScaler()),\n",
    "        (\"random_forest\", rf_clf),\n",
    "    ]\n",
    ")\n",
    "pipeline.fit(smiles_train, y_train)\n",
    "\n",
    "y_pred = pipeline.predict_proba(smiles_test)[:, 1]\n",
    "\n",
    "auroc = multioutput_auroc_score(y_test, y_pred)\n",
    "print(f\"Feature union AUROC for Random Forest: {auroc:.2%}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9352c7bdaf124c82",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "A major breakthrough in pretraining graph neural networks (GNNs) in chemistry was the paper [\"Strategies for Pre-training Graph Neural Networks\" W. Hu et al.](https://arxiv.org/abs/1905.12265). Their best AUROC on BACE was 84.5%. Therefore, we managed to outperform a complex, pretrained GNN using simple molecular fingerprints without any pretraining or sophisticated methods."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d47dcdc-fddd-49d8-bcf3-dfa4e590a2f3",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-10-26T13:30:28.813123200Z",
     "start_time": "2024-10-26T13:30:28.767117700Z"
    },
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## Bonus exercise\n",
    "\n",
    "Peptides are small proteins, with typically up to a few dozen aminoacids. While they are frequently analyzed as text sequences, they are, are all chemistry and biology, molecules built from atoms. Peptide function prediction is an important task in bioinformatics and peptide-based drugs.\n",
    "\n",
    "One dataset from this domain comes from Long Range Graph Benchmark (LRGB), which proposed a 10-task peptide function prediction dataset with over 15k peptide molecules. For details, see [\"Long Range Graph Benchmark\" V. Dwivedi et al.](https://arxiv.org/abs/2206.08164). Perform molecule classification on this dataset, using predetermined benchmark splits. It uses AUPRC (Average Precision, AP) as a metric.\n",
    "\n",
    "1. Load dataset and splits.\n",
    "2. Parse SMILES strings as molecules.\n",
    "3. Split molecules and labels. You can ignore the validation set.\n",
    "4. Calculate fingerprints, e.g. ECFP, Topological Torsion, EState.\n",
    "5. Train Random Forest classifier. Since the dataset is quite large, use 500 trees, instead of default 100.\n",
    "6. Calculate multioutput AUPRC score of the model.\n",
    "7. Compare results with paper.\n",
    "\n",
    "Use [scikit-fingerprints documentation](https://scikit-fingerprints.github.io/scikit-fingerprints/) as necessary. If you want to expand this further, you can perform hyperparameter tuning on the predetermined validation split."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a152e04d9c6f14e",
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from skfp.datasets.lrgb import load_peptides_func, load_lrgb_mol_splits\n",
    "from skfp.preprocessing import MolFromSmilesTransformer\n",
    "\n",
    "\n",
    "smiles_list, y = load_peptides_func()\n",
    "train_idxs, valid_idxs, test_idxs = load_lrgb_mol_splits(\"Peptides-func\")\n",
    "\n",
    "mol_from_smiles = MolFromSmilesTransformer(n_jobs=-1)\n",
    "mols = mol_from_smiles.transform(smiles_list)\n",
    "\n",
    "mols = np.array(mols)\n",
    "\n",
    "mols_train = mols[train_idxs]\n",
    "mols_test = mols[test_idxs]\n",
    "\n",
    "y_train = y[train_idxs]\n",
    "y_test = y[test_idxs]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58e516dc-7368-4a14-8596-6c34fab9cb5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from skfp.fingerprints import ECFPFingerprint, EStateFingerprint, TopologicalTorsionFingerprint\n",
    "from skfp.metrics import extract_pos_proba, multioutput_auprc_score\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "\n",
    "\n",
    "for fp_name, fp_cls in [\n",
    "    (\"ECFP\", ECFPFingerprint),\n",
    "    (\"EState\", EStateFingerprint),\n",
    "    (\"TopologicalTorsion\", TopologicalTorsionFingerprint),\n",
    "]:\n",
    "    fp = fp_cls(n_jobs=-1)\n",
    "    X_train = fp.transform(mols_train)\n",
    "    X_test = fp.transform(mols_test)\n",
    "\n",
    "    clf = RandomForestClassifier(n_estimators=500, n_jobs=-1, random_state=0)\n",
    "    clf.fit(X_train, y_train)\n",
    "\n",
    "    y_pred = clf.predict_proba(X_test)\n",
    "    y_pred = extract_pos_proba(y_pred)\n",
    "\n",
    "    auprc = multioutput_auprc_score(y_test, y_pred)\n",
    "    print(f\"{fp_name} AUPRC: {auprc:.2%}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d9773ce-2991-4e52-8b80-1d171b13acf4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}