lab1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lab 1: A Ligand-based Virtual Screening Pipeline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Review the lab material and go through the entire notebook. \n",
    "The lab contains 5 exercises for you to solve. The entire lab is worth 2.5% of your final grade and each exercise is worth 0.4% of your final grade. Going through the full notebook is worth 0.5% of your final grade.\n",
    "Any extra credit or bonus exercises are worth an additional 0.4%.\n",
    "\n",
    "Labs are due within 48 hours and can be submitted by email to the course instructor and TA."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Motivation of Virtual Screening\n",
    "\n",
    "Virtual screening (VS) is a computational framework to predict and prioritize compounds as drug candidates by simulating their interactions with a target receptor or predicting other properties relevant to drug discovery, thereby reducing the number of compounds that need to be experimentally tested.\n",
    "\n",
    "Virtual screening is much faster and cheaper than experimental screening; we can test a magnitude of $10^9$ to $10^{12}$ compounds per day (equivalent to at least ~3 years of experimental screening per day of virtual screening). \n",
    "\n",
    "Simulation of molecular dynamics or molecular docking are examples of VS approaches based on first principles of physical processes, such as intermolecular force fields, to discern how the compound or ligand might move around in physical space or reorient its structure to bind with the target."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading a Virtual Screening Library"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the following code blocks to set up the data for the lab. These will take some time to run so read on in the interim!\n",
    "\n",
    "We start by setting up installs according to the `requirements.txt` file and importing required modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install -r requirements.txt --quiet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import subprocess\n",
    "import time\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from rdkit import Chem\n",
    "from rdkit.Chem import (\n",
    "    AllChem,\n",
    "    Descriptors,\n",
    "    MolFromSmiles,\n",
    "    PandasTools,\n",
    ")\n",
    "\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)\n",
    "\n",
    "colors = [\"#A20025\", \"#6C8EBF\"]\n",
    "sns.set_palette(sns.color_palette(colors))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit.Chem import Draw\n",
    "d2d = Draw.MolDraw2DSVG(-1, -1)\n",
    "dopts = d2d.drawOptions()\n",
    "dopts.useBWAtomPalette()\n",
    "dopts.setHighlightColour((.635,.0,.145,.4))\n",
    "dopts.baseFontSize = 1.0\n",
    "dopts.additionalAtomLabelPadding = 0.15\n",
    "dopts.dotsPerAngstrom = 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit.Chem import PandasTools\n",
    "\n",
    "specs = PandasTools.LoadSDF(\n",
    "  \"data/Specs.sdf\", smilesName='smiles', molColName=None\n",
    ")[[\"PUBCHEM_SUBSTANCE_ID\", \"smiles\"]]\n",
    "\n",
    "# Note: This may take some time\n",
    "PandasTools.AddMoleculeColumnToFrame(specs, 'smiles', 'mol')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "specs.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Virtual Screening Taxonomy\n",
    "\n",
    "Virtual screening pipelines are often hierarchical, proceeding with a broad, large number of compounds through simpler, less restrictive, and computationally cheaper processes to whittle down the quantity of compounds for more complex processes. The simplest filters include compound filters for drug-likeness and known structural alerts. Since lead optimization results in an increase in molecular complexity, filters that assess lead-likeliness or desired ADMET properties may also be applied to discourage compounds that would necessitate too many changes during lead optimization.\n",
    "\n",
    "After use of chemical compound filters, virtual screening methods are broadly categorized depending on the information we have available about the target and its 3D structure. If we have the target’s 3D structure available, we can perform structure-based methods such as protein-ligand docking. By “3D structures,” we mean that previous experiments (using techniques such as X-ray crystallography or nuclear magnetic resonance spectroscopy) have determined the target’s structure.\n",
    "\n",
    "If a defined target structure is unavailable, either from a database such as PDB or from our own experiments, then we can still take advantage of ligand-based methods. If we have just one active compound, we can perform a similarity search: screen the active compound against a library of other compounds to identify which compounds are like the active compound. Intuitively, compounds with similar structures or properties to a known active compound have a higher chance of being active. If we have several known, active compounds, then we can also conduct a ligand-based pharmacophore search (in addition to a similarity search). In pharmacophore searching, we take our list of active compounds, generate possible configurations of their 3D shape, and try to align these possible configurations to find common structural features shared by the different active compounds. We can score structural features with high commonality and conduct a search for which of our unknown compounds have high-scoring features.\n",
    "\n",
    "If we have many known, active compounds as well as known, inactive compounds, then we can enter the realm of machine learning (ML) methods. In ligand-based machine learning, our goal is to train a model that can identify a correlation between the structural features of the compounds in our dataset and how they relate to the compound’s activity. If our model can discriminate between known active and inactive compounds, we hope it’s performance will generalize to new compounds that weren’t in the training dataset. Machine learning models for predicting compound activity are sometimes referred to as quantitative structure-activity relationship (QSAR) models (though not all QSAR models are ML models). Quantitative structure-property relationship (QSPR) is the more general term that refers to quantifying the relationship between structure and any molecular property.\n",
    "\n",
    "It is also common to use multiple virtual screening methods, one after the other, of increasing complexity. This approach is referred to as hierarchical virtual screening. Each step filters out more structures of no further interest, until at the end of the process a series of candidate structures remain. For example, we can use rule-based compound filters to rapidly screen out the most egregious candidates, followed by ligand-based virtual screening as a more nuanced detector of active compounds. Remaining compounds can then be subject to more rigorous and intensive structure-based methods, e.g., protein-ligand docking or molecular dynamics simulations, with one or more targets of interest."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![UN01](./figures/L01_UN01_Flynn.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've begun discussing data structures and algorithms to represent molecules, which are also handy for comparing molecules to each other (e.g., similarity) and sorting those molecules based on their similarity to a query structure. Assuming we have one active query compound, let's extend our knowledge to screen the query compound against a library of other compounds to identify which compounds are like the query."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scenario: Hit Identification of Antimalarial Compounds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We’ll show how to use ligand-based virtual screening to evaluate compounds for antimalarial potential based on the similarity of their structures to known antimalarial compounds. Malaria is a life-threatening parasitic disease that affects millions of people worldwide, primarily in tropical and subtropical regions. It is a major global health problem, causing a significant number of deaths, particularly among young children and pregnant women. Developing effective antimalarial drugs can directly contribute to reducing malaria’s global health burden and negative economic impact. \n",
    "Antimalarial drug research often focuses on various protein targets within Plasmodium falciparum, which is the most common vector of malaria. These protein targets are critical to its survival, replication, and invasion of host cells. An ominous threat to controlling malaria is P. falciparum’s ability to adapt and develop resistance to antimalarial drugs. Drugs that were once effective in treating malaria can become ineffective over time, and new therapeutics are repeatedly necessary.\n",
    "\n",
    "We will initiate a virtual screening campaign to identify antimalarial drug candidates, evoking recent successes for inspiration. We will begin with a sizable quantity of novel molecules and filter out substances that lack drug-like qualities or have detrimental properties, such as toxicity risk. To initiate the screen, we will start with a diverse range of >212K small molecules made available by PubChem from the SPECS repository, casting a wide net over chemical space. A key feature of SPECS is commercial availability – desirable compounds for downstream experiments are synthesizable and conveniently available for purchase. There are many other commercial vendors offering millions of compounds that meet a variety of characteristics and are purchasable, e.g., ChemBridge, eMolecules, Enamine, Life Chemicals, and Maybridge. \n",
    "\n",
    "We will then use a complementary resource, the Malaria Box compounds, containing 400 known actives against P. falciparum. To identify promising hit compounds, we will conduct a similarity search to cross-reference our diverse chemical library against the Malaria Box. In later chapters, we will consider how to optimize and develop discovered hits into potential antimalarial drugs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Strategy: Similarity Seaching & Sorting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar compounds have similar properties. If we have at least one compound that has a known activity or property of interest, we can use it as a reference point to search for similar compounds that might produce a similar effect. Naively, we could calculate similarity between the reference compound and each compound in our library, order the library compounds by highest similarity, and retain the most similar compounds for experimental validation. \n",
    "\n",
    "However, there are three key factors that add complexity to our naïve similarity searching strategy.\n",
    "1. Representation: How do we represent molecules as meaningful features to a computer? There are different ways we can represent a molecule in terms of features that describe its properties or structure.\n",
    "2. Similarity Metric: How do we compare similarity between two compounds? There are different metrics that quantify similarity between molecules as a proxy for similar biological function.\n",
    "3. Search Strategy: If there is one reference, we may only care about the most similar molecule(s). But if there are multiple references, how should we prioritize a molecule based on its aggregate similarity to all available reference compounds? What if we want a search strategy that can accommodate new molecules that we might consider later? Furthermore, molecules can be similar in different ways. How do we factor this into our search? Which features are meaningful for similarity comparison?\n",
    "\n",
    "Each factor has multiple available options with advantages or disadvantages, and they are linked together. In general, our similarity metric measures the extent of shared features between molecule pairs. Choice of feature representation affects results of using different similarity metrics and search strategies. As an example, the following figure exemplifies the nuances in structure-property relationships that make similarity searching a challenging task."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![UN02](./figures/L01_UN02_Flynn.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Structure Data Files (SDFs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We loaded the Specs data set by ingesting it as a Structure Data File (SDF) into a Pandas DataFrame. No need to have familiarity with Pandas, it is just a tool we can use to store the molecules in our dataset.\n",
    "\n",
    "SDFs are one way we can store multiple molecular structures (both 2D and 3D) along with associated property information. SDFs are derived from the concept of a connection table. The simplest connection table defines enough information to generate a molecule as a graph, where nodes are atoms and edges are bonds. This only requires an atom table and a bond table.\n",
    "\n",
    "An atom table indexes each atom. A bond table designates the two atoms involved in each bond and the type of bond. Bond orders of 1, 2, and 3 correspond to single bonds, double bonds, and triple bonds, respectively. Hydrogens may be implicitly or explicitly defined. The former is not an issue as algorithms exist that can determine the number of hydrogens. Also note that connection tables are not unique – we could derive a different, but equivalent connection table by swapping the indices of two atoms and their bonds.\n",
    "To express meaningful chemical information cohesively and explicitly, SDFs extend connection tables to incorporate 2D or 3D spatial information and information associated with individual atoms or bonds and the entire structure. The most common SDF format is based on the MOL file format developed by Molecular Design Limited (MDL), of which the V2000 and V3000 formats are often used.\n",
    "\n",
    "Information for each SDF entry is compartmentalized into blocks. The bottom-half of the following figure shows an example entry for benzoic acid. The header block’s two lines describe the molecule’s name or formula (C7H6O2), the program used to make it (ChemDraw), the date and time it was made (August 15th, 2023), and if 2D or 3D coordinates are given. The count block tells us there are 9 atoms, 9 bonds, no chirality, and that the file format is V2000. The atom block contains X, Y, Z coordinates of each atom and the atom symbol. The additional columns may encode property information, but a separate properties block is usually used instead. The bonds block designates two atoms by their index (first and second columns) and their bond order (third column) and stereochemistry (fourth column)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![UN03](./figures/L01_UN03_Flynn.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's inspect one of the molecules and its type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mol = MolFromSmiles(specs.smiles[0])\n",
    "mol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The following prints out a cleaner version of the above\n",
    "img = Draw.MolsToGridImage(\n",
    "  mols=[mol], molsPerRow=1, useSVG=True, drawOptions=dopts,\n",
    ")\n",
    "img"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(mol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Molecule Sanitization\n",
    "\n",
    "As the Pandas DataFrame was being processed, you likely noticed a few errors with the form “ERROR: Could not sanitize molecule ending on line 6498564” followed by an explanation of why the SMILES produced the error. These errors are related to sanitization, a preprocessing step that validates whether the SMILES are chemically reasonable before converting it into a Mol. \n",
    "\n",
    "Sanitization enforces data quality and reliability to avoid incorrect conclusions, wasted resources, and flawed research outcomes. What constitutes a molecule as valid or not depends on the use case. In virtual screening, we may remove compounds with highly reactive motifs that may manifest as off-target effects or we may remove compounds with bad elements in the context of medicinal chemistry, e.g., compounds with a metal atom. \n",
    "\n",
    "At a minimum, we wouldn’t want to retain compounds that aren’t chemically sensible or break fundamental constraints such as valence rules. For most data sets we work with in part one, we will assume compounds are chemically valid if they pass RDKit’s default validation and sanitization. Optionally, RDKit sanitization can be disabled, e.g., mol = Chem.MolFromSmiles(smiles, sanitize=False).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Structuring the Pipeline\n",
    "\n",
    "Process:\n",
    "- Apply compound filters (new topic!)\n",
    "  - Compute molecular descriptors\n",
    "  - Property-based Filters\n",
    "  - Structure-based Filters\n",
    "- Fingerprinting our Library\n",
    "- Similarity Searching\n",
    "- Sorting (Ranking)\n",
    "\n",
    "## Molecular Descriptors\n",
    "\n",
    "We want to explore and filter out molecules based on quantifiable properties. We currently have the SMILES column available as a feature. But SMILES are a text string, not a number. The SMILES need to be processed in a way that allows a filter to evaluate it or a model to learn from it. Molecules are complex and how we process them affects a model’s ability to learn from them. There has been a lot of research on how to quantify or featurize molecules, and we will review popular approaches. First, we will compute molecular descriptors that can be used for compound filtering. Later in this chapter, we will compute fingerprint descriptors that can be used for similarity searching.\n",
    "\n",
    "Molecular descriptors are useful for interpreting a molecule’s properties and for encoding a molecule’s chemical information, so that we can develop a model to predict the properties of other molecules. There are thousands of descriptors that quantify different molecular features and that can be calculated using many different software tools. There are additional, experimentally derived descriptors that can also be accessed, such as databases of empirical molecular constants. Naturally, there is no shortage of documentation on descriptors, including the handbook on molecular descriptors. The below table summarizes the taxonomy of molecular descriptors. RDKit can compute 211 descriptors; let’s dive-in and calculate descriptors on our dataset.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![UN04](./figures/L01_UN04_Flynn.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "RDKit's `Descriptors` library provides implementations for 211 descriptors. We print out the first 3, which provides both the name of the descriptor and its function call."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit.Chem import Descriptors\n",
    "print(len(Descriptors._descList))\n",
    "print(Descriptors._descList[:3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code block shows how we can use RDKit to compute 4 descriptors that describe each molecule’s molecular weight, number of hydrogen bond acceptors (HBAs), number of hydrogen bond donors (HBDs), and logP. \n",
    "\n",
    "These descriptors will be useful to define a filter that removes compounds that are not \"drug-like? What does it mean to be \"dru-like\"? Well, given a large library of compounds, we might want some simple rules to characterize them as either drug-like or not drug-like. We could then use our rules for defining drug likeness to filter out compounds that are not drug-like, leaving us to focus on a smaller set of drug-like compounds.\n",
    "There are multiple guidelines for classifying compounds by their drug-likeness. Perhaps the most common and simple among them is Lipinski’s Rules or Rule of Five, which states that a compound is not drug-like if it violates more than one of the following criteria:\n",
    "\n",
    "1.\tMolecular weight less than 500 daltons\n",
    "2.\tLipophilicity (calculated logP) of less than 5\n",
    "3.\tNo more than 10 hydrogen bond acceptors\n",
    "4.\tNo more than 5 hydrogen bond donors\n",
    "\n",
    "Note that all the rules specify quantity limits that are multiples of five, hence Lipinski’s Rule of Five. Though beneficial for their simplicity, these rules make a crude prediction of a compound’s potential as an oral drug based on properties that affect its bioavailability. Depending on date of analysis, approximately 16% to 20% of orally administered drugs on the market violate this guideline \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RDKIT_DESCRIPTORS = { desc : func for desc, func in Descriptors._descList }\n",
    "RO5_PROPS = ['ExactMolWt', 'NumHAcceptors', 'NumHDonors', 'MolLogP']\n",
    "\n",
    "def compute_descriptors(mol, func, missing_val=None):\n",
    "  try:\n",
    "    return func(mol)\n",
    "  except:\n",
    "    return missing_val\n",
    "\n",
    "for desc in RO5_PROPS:\n",
    "  specs[desc] = specs[\"mol\"].apply(\n",
    "    lambda x: compute_descriptors(x, RDKIT_DESCRIPTORS[desc])\n",
    "  )\n",
    "\n",
    "specs = specs.dropna(subset=RO5_PROPS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the `head()` and `describe()` functions to explore what the top 5 rows of the data set look like and to summarize basic statistics of the data set's numerical properties, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "specs.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "specs.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compound Filters\n",
    "\n",
    "Filtering out compounds is a low-cost, scalable maneuver to cut the fat from a large screening library and save downstream effort and resources. We can think of compound filters across two dimensions: property-based filters that consider the values of the descriptors we calculated in the previous section and substructure-based filters that consider the type and frequency of the substructures that we can decompose our molecules into. \n",
    "\n",
    "### Property-based Filters\n",
    "\n",
    "A convenient property-based filter is Lipinski’s Rule of Five (Ro5) as a starting point for evaluating drug-likeness. Assessing a compound’s adherence to Ro5’s four simple criteria measures prospect of oral bioavailability. Molecules that violate Ro5 criteria are likely to constitute poor solubility and permeability characteristics that hinder the effectiveness of oral drugs. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lipinski_filter(row):\n",
    "    \"\"\"\n",
    "    Apply Lipinski's Rule of Five to filter molecules.\n",
    "    Returns True if the molecule violates no more than one rule.\n",
    "    \"\"\"\n",
    "    violations = 0\n",
    "    # Check each rule\n",
    "    if row[3] > 500: violations += 1\n",
    "    if row[6] > 5: violations += 1\n",
    "    if row[5] > 5: violations += 1\n",
    "    if row[4] > 10: violations += 1\n",
    "\n",
    "    # Return True if no more than one rule is violated\n",
    "    return violations <= 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "specs[\"ro5_compliant\"] = specs.apply(lipinski_filter, axis=1)\n",
    "specs_ro5_compliant = specs[specs[\"ro5_compliant\"]]\n",
    "specs_ro5_violated = specs[~specs[\"ro5_compliant\"]]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After applying the R05 filter, we've reduced the compound library from 212,670 compounds to 199,481 \"drug-like\" compounds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Compound library size pre-RO5 filter: {len(specs)}\")\n",
    "print(f\"Compound library size post-RO5 filter: {len(specs_ro5_compliant)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice also that the distribution of the descriptors has changed. For example, high molecular weight compounds have been filtered out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "specs_ro5_compliant.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![UN05](./figures/L01_UN05_Flynn.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The other 3 descriptors are within closer numerical ranges to each other, so we visualize their changes in one lollipop plot comparing the average values for number of HBAs, logP, and number of HBDs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ro5_avgs_df = pd.DataFrame(data = {\n",
    "    'property' : ['NumHAcceptors', 'NumHDonors', 'MolLogP'],\n",
    "    'RO5 Compliant Mols' : [specs_ro5_compliant['NumHAcceptors'].mean(), specs_ro5_compliant['NumHDonors'].mean(), specs_ro5_compliant['MolLogP'].mean()],\n",
    "    'RO5 Violated Mols' : [specs_ro5_violated['NumHAcceptors'].mean(), specs_ro5_violated['NumHDonors'].mean(), specs_ro5_violated['MolLogP'].mean()]\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ro5_avgs_df['rel_change'] = (ro5_avgs_df['RO5 Compliant Mols'] - ro5_avgs_df['RO5 Violated Mols']) / ro5_avgs_df['RO5 Violated Mols']\n",
    "ordered_df = ro5_avgs_df.sort_values(by='RO5 Compliant Mols')\n",
    "my_range = range(1, len(ordered_df.index) + 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![UN06](./figures/L01_UN06_Flynn.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Structure-based Filters\n",
    "\n",
    "Substructure-based filters are another tool for early elimination of compounds that are unlikely to meet desired criteria. Displeasing substructures interfere with experimental assays or have undesirable properties (e.g., are toxic).\n",
    "\n",
    "#### Pan Assay Interference Compounds (PAINS)\n",
    "\n",
    "Promiscuous bioactive compounds, or PAINS compounds, appear active in many high-throughput screening experiments against a broad, diverse range of targets. PAINS are often false positive hits and mislead results in biological assays, producing artifacts in the data. These “frequent hitters” exhibit a propensity to bind to multiple targets due to nonspecific binding or interference with the specific assay technology readouts. PAINS filters are composed of substructures that are associated with PAINS. Conveniently, RDKit maintains pre-defined PAINS filters we can easily use to filter our compound library.\n",
    "\n",
    "#### Structural Alerts for Toxicity\n",
    "\n",
    "Broadly, structural alerts refer to any filter based on a molecule’s structural composition. In this context, we will assemble a set of structural alerts that have a known association with toxic endpoints. For example, thiophene, epoxides, and acid anhydrides are common structural alerts used to identify risk of hepatotoxicity, mutagenicity, or skin sensitization, respectively. If a compound contains several structural alerts greater than a user-defined threshold, it is filtered out.\n",
    "\n",
    "Cross-referencing compounds to a list of structural alerts is a widely used, simple method for identifying problematic compounds. However, structural alerts have limited prognostic utility – they are not sufficient to declare a molecule as toxic and they can misclassify toxic molecules as safe. Since structural alerts are determined retrospectively, they have limited predictive power for new or understudied substructures. \n",
    "\n",
    "![UN07](./figures/L01_UN07_Flynn.svg)\n",
    "\n",
    "In addition to PAINS filters, we’ll also use BRENK filters, which were constructed to remove undesirable substructures that pose toxicity risks within the context of drug discovery for neglected diseases "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit.Chem import FilterCatalog\n",
    "filter_params = FilterCatalog.FilterCatalogParams()\n",
    "filter_params.AddCatalog(filter_params.FilterCatalogs.PAINS)\n",
    "filter_params.AddCatalog(filter_params.FilterCatalogs.BRENK)\n",
    "catalog = FilterCatalog.FilterCatalog(filter_params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "specs_ro5_compliant[\"PAINS_BRENK_compliant\"] = specs_ro5_compliant['mol'].apply(catalog.HasMatch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "specs_ro5_pains_brenk_compliant = specs_ro5_compliant[specs_ro5_compliant[\"PAINS_BRENK_compliant\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After successive application of the RO5, PAINS, and BRENK filters, we're down to 93,915 compounds. That is less than half of our starting compound library size!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"SPECS compounds compliant with RO5, PAINS, and BRENK: {specs_ro5_compliant['PAINS_BRENK_compliant'].sum()}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1: Compound Filters\n",
    "\n",
    "In the code walkthrough, we used BRENK and PAINS filters from the RDKit catalog. Pick another filter from the RDKit catalog and apply it to the compound library (`specs_ro5_compliant`) instead of the BRENK and PAINS filters. What type of structures or properties does the filter you chose assess? How many compounds remain after applying the filter? Was it a more liberal or conservative filter than the ones we used?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Student Solution to Exercise 1\n",
    "\n",
    "Provide your solution to the above exercise in this cell and/or immediately following cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: Student to provide solution below"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Rapid Elimination of Swill (REOS)\n",
    "\n",
    "REOS quickly filters out low quality compounds, i.e., “swill”, that are unlikely to be viable drug candidates. REOS improves screening efficiency by avoiding pursuit of unproductive compounds with high likelihood of causing assay artifacts or false positives. RDKit does not have a pre-defined filter for the for REOS, as REOS is more of a concept than a specific ruleset. For example, PAINS is a specific ruleset that embodies REOS’ goal. But that does not block us from defining and using a custom substructure filter with RDKit’s HasSubstructMatch() function.\n",
    "\n",
    "As an example, let’s load a set of publicly available structure alerts that define the Glaxo Wellcome hard filters. This filter is conveniently encoded as SMARTS. SMARTS are an extension of SMILES that represent specific substructural patterns. We will use SMARTS to search for substructures that might be within each molecule in our library, like regular expressions adapted to chemical compounds. We recommend Daylight’s reference page for a more information: https://www.daylight.com/dayhtml/doc/theory/theory.smarts.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit.Chem import Draw, MolFromSmarts\n",
    "glaxo_alerts = pd.read_csv(\"data/glaxo_structural_alerts.csv\")\n",
    "glaxo_alerts[\"ROMol\"] = glaxo_alerts.smarts.apply(MolFromSmarts)\n",
    "print(f\"Number of Glaxo Wellcome structural alerts: {len(glaxo_alerts)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We visualize examples of 6 substructures within the set of Glaxo Wellcome hard filters. Compounds in our library that have these substructures will be filtered out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = Draw.MolsToGridImage(\n",
    "  mols=glaxo_alerts[\"ROMol\"].iloc[2:8].tolist(), molsPerRow=3, legends=glaxo_alerts[\"description\"].iloc[2:8].tolist(), useSVG=True, drawOptions=dopts,\n",
    ")\n",
    "img"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Applying the Glaxo Wellcome hard filter sifts out another chunk of molecules, leaving us with 61,372 compounds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glaxo_sa_matches = []\n",
    "def glaxo_filter(mol, alerts):\n",
    "  match_ = False\n",
    "  for _, alert in alerts.iterrows():\n",
    "    if mol.HasSubstructMatch(alert.ROMol):\n",
    "      glaxo_sa_matches.append({\n",
    "        \"mol\": mol,\n",
    "        \"alert\": alert.ROMol,\n",
    "        \"description\": alert.description,\n",
    "      })\n",
    "      match_ = True\n",
    "  return match_\n",
    "\n",
    "specs_ro5_pains_brenk_compliant[\"GLAXO_compliant\"] = ~specs_ro5_pains_brenk_compliant['mol'].apply(glaxo_filter, alerts=glaxo_alerts)\n",
    "specs_filtered = specs_ro5_pains_brenk_compliant[specs_ro5_pains_brenk_compliant[\"GLAXO_compliant\"]]\n",
    "glaxo_sa_matches = pd.DataFrame(glaxo_sa_matches)\n",
    "print(f\"Compound library size after filtering: {specs_ro5_pains_brenk_compliant['GLAXO_compliant'].sum()}\")\n",
    "glaxo_sa_matches[\"description\"].value_counts()[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We show a few examples of compound matches to the Glaxo Wellcome filter, with the match substructures highlighted in light red."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d2d_hl = Draw.MolDraw2DSVG(-1, -1)\n",
    "dopts_hl = d2d_hl.drawOptions()\n",
    "dopts_hl.useBWAtomPalette()\n",
    "dopts_hl.setHighlightColour((.635,.0,.145,.5))\n",
    "dopts_hl.baseFontSize = 1.0\n",
    "dopts_hl.additionalAtomLabelPadding = 0.15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "highlights = [mol.GetSubstructMatch(alert) for mol, alert in zip(glaxo_sa_matches.mol, glaxo_sa_matches.alert)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = Draw.MolsToGridImage(\n",
    "  glaxo_sa_matches.mol.iloc[2:8].tolist(),\n",
    "  highlightAtomLists=highlights[2:8],\n",
    "  molsPerRow=3,\n",
    "  legends=glaxo_sa_matches.description.iloc[2:8].tolist(),\n",
    "  drawOptions=dopts_hl,\n",
    "  useSVG=True\n",
    ")\n",
    "img"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With just a few simple filters, we’ve reduced our starting library to a third of its original size without losing any molecule’s worth further investigation. Remaining downstream processes will extract greater computational cost, so we’ve reduced noise and improved efficiency of the overall pipeline we are constructing. Let’s move on to the next step, which entails numerically representing our molecules as fingerprints, which we can use to conduct a similarity search and identify promising antimalarial candidates."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fingerprints: Representing Molecules as Numbers\n",
    "\n",
    "### Fingerprint our Library"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For our task, we’ll stick with Morgan (circular) fingerprints to generate fingerprint descriptors. Note that, in practice, it is common to try multiple types of fingerprints to assess performance for a given task. To aid in reproducibility, you should always state what fingerprints you used to featurize a molecule and the parameters. We will use a radius of 2 and a fingerprint of length 2048.\n",
    "\n",
    "Why do we fingerprint our library with a radius of 2 and not 4? We previously described one benefit of radius 4 fingerprints is their ability to capture a larger molecular context, providing a more global perspective on the molecular structure. Features involving atoms that are further apart from each other are more likely to be included. Inversely, a smaller radius of 2 highlights the immediate neighborhood of each atom and is more likely to capture specific functional groups, local arrangements of atoms, or other fine-grained details of molecular structure. Calculating Morgan fingerprints with a smaller radius is computationally cheaper, which is pedagogically beneficial since it is easier to modify the example code and produce new results. In practice, we might try and compare both."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit import DataStructs\n",
    "from rdkit.Chem import AllChem\n",
    "\n",
    "def compute_fingerprint(mol, r, nBits):\n",
    "    fp = AllChem.GetMorganFingerprintAsBitVect(mol, r, nBits=nBits)\n",
    "    return fp\n",
    "\n",
    "specs_filtered[\"morgan_fp_r2_b2048\"] = specs_filtered[\"mol\"].apply(compute_fingerprint, r=2, nBits=2048)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code blocks show how we can decompose a molecule and visualize the bits that are set for its corresponding fingerprint. In this example, we only show 5 bits, but the code is easily altered to show all bits for any molecule (this particular molecule has a total of 47 bits set). Notice how none of the fragments cover a path of more than 2-steps away from their central atom, which is highlighted by a blue circle. Yellow circles indicate aromatic atoms; grey circles indicate aliphatic ring atoms. Atoms or bonds that are light grey indicate structures that influence connectivity but are not directly part of the fingerprint. Asterisks indicate extensions beyond the fingerprint itself and can be interpreted as wildcards."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_fragment_from_bit(mol, bit_number):\n",
    "  \"\"\" Given an rdkit mol, draws the local fragment that corresponds to the set bit of ecfp featurization.\n",
    "\n",
    "  If the bit is not set, will throw an error.\n",
    "  \"\"\"\n",
    "  bi = {}\n",
    "  fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius=2, bitInfo=bi)\n",
    "  try:\n",
    "    svg = Draw.DrawMorganBit(mol, bit_number, bi, useSVG=True)\n",
    "  except:\n",
    "    raise ValueError(f\"Featurization of mol doesn't have bit {bit_number} set\")\n",
    "  return svg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "example_mol = specs_filtered.mol.iloc[0]\n",
    "img = Draw.MolsToGridImage(\n",
    "  mols=[example_mol], molsPerRow=1, useSVG=True, drawOptions=dopts,\n",
    ")\n",
    "img"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "example_fp = specs_filtered.morgan_fp_r2_b2048.iloc[0]\n",
    "fp_hit_indices = [idx for idx, bit in enumerate(example_fp) if bit][:5]\n",
    "example_bits = [draw_fragment_from_bit(example_mol, hit_idx) for hit_idx in fp_hit_indices]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(*example_bits)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We package this into a nice figure for ease of viewing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![UN08](./figures/L01_UN08_Flynn.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Similarity Searching\n",
    "\n",
    "### Defining \"Similarity\"\n",
    "\n",
    "If we take two molecules, we can think of their fingerprints as representing the coordinates of two different data points in a fingerprint feature space. Tanimoto similarity and Dice similarity are two similarity coefficients commonly used to quantify the similarity between sets or binary data, e.g., fingerprints.\n",
    "\n",
    "Tanimoto similarity (also known as Jaccard similarity) is a measure of similarity between two sets, A and B, and is calculated as the ratio of the size of their intersection to the size of their union. In the context of binary vectors, Tanimoto similarity is calculated as the number of common non-zero elements divided by the total number of non-zero elements in both vectors.\n",
    "\n",
    "Dice similarity is another measure of similarity between two sets, A and B. It's calculated as twice the size of their intersection divided by the sum of the sizes of the individual sets. For binary vectors, Dice similarity measures the proportion of common non-zero elements relative to the total number of non-zero elements summed across both vectors. \n",
    "Mathematically, Tanimoto similarity (T) and Dice similarity (D) are illustrated in the figure below.\n",
    "\n",
    "Both Tanimoto similarity and Dice similarity range from 0 to 1, where 0 indicates no similarity and 1 indicates complete similarity between the sets or vectors being compared. Note that a similarity of 1 does not indicate that the two molecules are the same. The Daylight manual provides a quick reference to other common similarity measures.\n",
    "\n",
    "![UN09](./figures/L01_UN09_Flynn.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we commence a similarity search across all 400 Malaria Box compounds, we familiarize ourselves with the behavior of Tanimoto and Dice similarities computed between our remaining compound library and a single Malaria Box compound as a query. The following figure shows the extremes – the 3 most similar molecules to the query and the 3 least similar molecules. We can clearly notice a difference in the structures for the 3 most similar molecules and the 3 least similar molecules, where the latter look nothing like the query. In this case, Tanimoto and Dice similarity metrics result in the same molecules at the extremes, though the magnitude of similarity differs. From these few examples, it seems that the Dice coefficient results in higher similarity scores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "malaria_box = pd.read_excel(\"data/MalariaBox400compoundsDec2014.xls\", usecols=[\"HEOS_COMPOUND_ID\", \"Smiles\"])\n",
    "PandasTools.AddMoleculeColumnToFrame(malaria_box, 'Smiles', 'mol')\n",
    "malaria_box[\"morgan_fp_r2_b2048\"] = malaria_box[\"mol\"].apply(compute_fingerprint, r=2, nBits=2048)\n",
    "\n",
    "query = malaria_box.morgan_fp_r2_b2048.iloc[236]\n",
    "mols = specs_filtered.morgan_fp_r2_b2048.tolist()\n",
    "\n",
    "specs_filtered[\"tanimoto_sim\"] = DataStructs.BulkTanimotoSimilarity(query, mols)\n",
    "specs_filtered[\"dice_sim\"] = DataStructs.BulkDiceSimilarity(query, mols)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![UN10](./figures/L01_UN10_Flynn.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2: Fingerprint Similarity Metrics\n",
    "\n",
    "We computed Tanimoto and Dice similarity coefficients. Pick another similarity coefficient mentioned in the referenced article from the Daylight manual: https://www.daylight.com/dayhtml/doc/theory/theory.finger.html (see section 6.3).\n",
    "\n",
    "How does your chosen similarity coefficient differ from Tanimoto and Dice similarity? Does it result in similarity scores that tend to be greater or smaller than Tanimoto similarity? You can modify the plotting code below to visualize the relative difference w.r.t. Tanimoto similarity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Student Solution to Exercise 2\n",
    "\n",
    "Provide your solution to the above exercise in this cell and/or immediately following cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: Student to provide solution below\n",
    "\n",
    "# Plotting code you can modify and re-use (swap out \"dice_sim\" for your chosen similarity metric)\n",
    "plt.figure()\n",
    "sns.pairplot(specs_filtered[[\"tanimoto_sim\", \"dice_sim\"]])\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Influence of similarity metric on scoring\n",
    "\n",
    "Choice of a similarity score threshold is dependent on the similarity measure. We can plot the distribution of both Tanimoto and Dice similarity scores against each other to check if this trend generalizes. The below figure shows the relationship is nearly linear and confirms that the Dice similarity score is greater than or equal to Tanimoto similarity score for each molecule in our library. If we were to choose a similarity threshold, e.g., 0.65, we may expect that more molecules will register as hits in our similarity search if we use the Dice similarity metric instead of the Tanimoto similarity metric. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![UN11](./figures/L01_UN11_Flynn.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Beyond dependence on the similarity metric used, similarity search results are also highly dependent on the chosen descriptors or features used to represent each molecule. These factors can cause the search to be too conservative, only returning compounds that are too similar to the reference compound, or too loose such that too many false positive hits are retained. Regarding the former, hits that are so similar that they share common substructures with our reference compounds may not be patentability.\n",
    "\n",
    "One rule-of-thumb is that a Tanimoto similarity of 0.85 corresponds to equivalent biological activity for drug-like molecules. Depending on context of the overall drug discovery project or the specific biological process being modeled, acceptable thresholds can be lower (e.g., 0.65) or higher (e.g., 0.9) and we recommend consulting related literature for a given problem domain to understand precedence set by prior work. In our case, we’ll count a hit as a molecule in our library that has at least one match against a query with a Dice similarity score of at least 0.65. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our similarity search (below) discovered 1,892 suitable molecules based on our similarity criteria. Commonly, we may be limited further by a pre-defined budget of how many molecules we can afford to screen in downstream processes. Let’s suppose we our limited by a budget of 1000 molecules and conduct a brute-force search to retain the top 1000 molecules by similarity score."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import defaultdict\n",
    "specs_filtered = specs_filtered.reset_index(drop=True)\n",
    "mols = specs_filtered.morgan_fp_r2_b2048.tolist()\n",
    "matches = defaultdict(int)\n",
    "\n",
    "for query in malaria_box.morgan_fp_r2_b2048:\n",
    " dice_sim = DataStructs.BulkDiceSimilarity(query, mols)\n",
    " matches.update({\n",
    "    idx : max(sim, matches[idx]) for idx, sim in enumerate(dice_sim) if sim >= 0.65\n",
    " })\n",
    "\n",
    "print(f\"Total hits: {len(matches)}\")\n",
    "print(f\"Hits per query: {len(matches) / malaria_box.shape[0]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "brute_force_matches = []\n",
    "budget = 1000\n",
    "for idx, sim in matches.items():\n",
    "    if len(brute_force_matches) <= budget:\n",
    "        brute_force_matches.append((sim, idx))\n",
    "    else:\n",
    "        min_tm_idx = len(brute_force_matches)\n",
    "        min_sim = sim\n",
    "        for tm_idx, (sim2, idx2) in enumerate(brute_force_matches):\n",
    "            if sim2 < min_sim:\n",
    "                min_sim = sim2\n",
    "                min_tm_idx = tm_idx\n",
    "        if min_tm_idx < len(brute_force_matches):\n",
    "            brute_force_matches[min_tm_idx] = (sim, idx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "top_matches = [idx for sim, idx in brute_force_matches]\n",
    "top_matches"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To save our top hits, we can run the following code blocks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "specs_hits_to_malaria_box = specs_filtered.filter(items=top_matches, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "specs_hits_to_malaria_box.to_csv(\"specs_hits_to_malaria_box.csv\", columns=[\"PUBCHEM_SUBSTANCE_ID\", \"smiles\"], index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 3: Being an Intrepid Programmer!\n",
    "\n",
    "The brute-force implementation is really slow when the number of entries in `matches` greatly exceeds the `budget`. It is also prohibitively memory-intensive if our `budget` is large.\n",
    "\n",
    "Many programmers have had to deal with similar problems, and have developed methods or libraries to deal with their problems. We can take advantage of this be exploring well-known libraries available in our programming language's ecosystem. To improve our implementation, we can use the `heapq` library: https://docs.python.org/3/library/heapq.html\n",
    "\n",
    "1. Summarize what the heapq library accomplishes. What are 3 methods available in the library? What are the performance considerations (e.g., Big Oh) of these methods? Take an example code segment from the documentation page and explain what it is doing. \n",
    "2. Use `heapq` as an alternative method to sample the top 1000 hits.\n",
    "3. Compare runtime complexity and memory complexity of the brute-force implementation and your `heapq` implementation.\n",
    "4. Are the top 1000 compounds of your implementation the same as the top 1000 compounds from the brute-force implementation?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Student Solution to Exercise 3\n",
    "\n",
    "Provide your solution to the above exercise in this cell and/or immediately following cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import heapq\n",
    "budget = 1000\n",
    "# TODO: Student to implement solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bonus Exercise: Quickselect\n",
    "\n",
    "There is an algorithm, Quickselect, that we can modify to achieve the optimal runtime complexity and space complexity. Explain the Quickselect algorithm, what runtime and space Big Oh complexity we can achieve, and then implement a function that uses Quickselect to sample the top 1000 hits."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Student Solution to Bonus Exercise\n",
    "\n",
    "Provide your solution to the above exercise in this cell and/or immediately following cells."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Retrospective\n",
    "\n",
    "We only needed at least one known active molecule for antimalarial activity to enable ultra-fast screening via similarity search. However, even with a liberal similarity threshold, we only retained 1892 molecules as hits against the Malaria Box compounds. A starting library of over 212K small molecules seems much smaller once we realize how few hits might result after application of property, substructure, and similarity filters.\n",
    "\n",
    "Lastly, chemical similarity searching alone is not sufficient – it does not guarantee equivalent biological effect. In particular, activity cliffs refer to minor changes between two molecular structures that are associated with large differences in function. We'll revisit this point later in the course. \n",
    "\n",
    "Note: The numbers in the figure below represent a ballpark approximation of the actual amount of remaining compounds from running this lab.\n",
    "\n",
    "![UN12](./figures/L01_UN12_Flynn.svg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 4: How could we improve?\n",
    "\n",
    "As we saw in lecture, many problems have overlapping solution sets. Consulting external resources gives us ideas on how to improve. Pick one of the following papers to read and provide a half-page summary on what lessons from the paper could be applicable to our problem.\n",
    "\n",
    "- Roy, K.K. (2017). Targeting the active sites of malarial proteases for antimalarial drug discovery: approaches, progress and challenges. International Journal of Antimicrobial Agents, 50 (3): 287-302\n",
    "- Shibeshi, M.A., Kifle, Z.D., & Atnafie, S.A. (2020). Antimalarial Drug Resistance and Novel Targets for Antimalarial Drug Discovery. Infection and drug resistance, 13, 4047–4060. https://doi.org/10.2147/IDR.S279433\n",
    "- Brenk, R., Schipani, A., James, D., Krasowski, A., Gilbert, I.H., Frearson, J., & Wyatt, P.G. (2008). Lessons learnt from assembling screening libraries for drug discovery for neglected diseases. Chem. Med. Chem., 3 (3): 435–444. https://doi.org/10.1002/cmdc.200700139\n",
    "- Walters, W.P., Stahl, M.T., Murcko, M.A. (1998). Virtual Screening -- An Overview. Drug Discov. Today. 3: 160–178"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Student Solution to Exercise 4\n",
    "\n",
    "Provide your solution to the above exercise in this cell and/or immediately following cells."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 5: Exploration\n",
    "\n",
    "Find another compound data set different from the Malaria Box compounds. If you conduct a similarity search of the Specs compound library (post compound filtering) against your new data set, do you get any hits? Do you get more or less hits (total and per query) compared to our reference against the Malaria Box compounds? Describe the data set you chose and its potential applications in two sentences or less."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Student Solution to Exercise 5\n",
    "\n",
    "Provide your solution to the above exercise in this cell and/or immediately following cells."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ml4dd",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}