02_rdkit_intro.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c20abcbf-1e4c-4952-b6c7-1e3c7ed2721a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# apply Jupyter notebook style\n",
    "from IPython.core.display import HTML\n",
    "\n",
    "from custom.styles import style_string\n",
    "\n",
    "HTML(style_string)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9020474",
   "metadata": {},
   "source": [
    "<div style=\"text-align:center;\">\n",
    "  <img src=\"custom/molssi_main_horizontal.png\" style=\"display: block; margin: 0 auto; max-height:200px;\">\n",
    "</div>\n",
    "\n",
    "\n",
    "# Introduction to RDKit\n",
    "\n",
    "<div class=\"overview admonition\"> \n",
    "<p class=\"admonition-title\">Overview</p>\n",
    "\n",
    "Questions:\n",
    "\n",
    "* What is RDKit?\n",
    "\n",
    "* How can I get molecule information in Python?\n",
    "\n",
    "* What is a molecular descriptor?\n",
    "\n",
    "Objectives:\n",
    "\n",
    "* Use RDKit to create molecules in Python\n",
    "\n",
    "</div>\n",
    "\n",
    "There are Python libraries that are made for working just with chemical data. One commonly used library in Python for data science (or cheminformatics) is called [RDKit](https://en.wikipedia.org/wiki/RDKit). RDKit is an open-source cheminformatics library, primarily developed in C++, and was created by [Dr. Greg Landrum](https://scholar.google.com/citations?user=xr9paY0AAAAJ&hl=en) in 2006. We will be using the Python interface to RDKit, though there are interfaces in other languages.\n",
    "\n",
    "The primary goal of RDKit is to facilitate the manipulation of chemical structures, including reading and writing molecular file formats, calculating molecular properties, and performing substructure searches. In addition, it offers a wide range of cheminformatics algorithms such as molecular fingerprint generation, similarity metrics calculation, and molecular descriptor computation. This notebook introduces RDKit basis, but we will see more of all of these topics in later notebooks.\n",
    "\n",
    "<div class=\"tip admonition\"> \n",
    "<p class=\"admonition-title\">Python Objects</p>\n",
    "\n",
    "Most of this functionality is achieved through the RDKit `mol` object. In Python, we use the word \"object\" to refer to a variable type with associated data and methods. \n",
    "One example of an object we have seen in notebooks is a list - we could also call it a \"list object\". An object has `attributes` (data) and `methods`. \n",
    "You access information about objects with the syntax\n",
    "```python\n",
    "object.data\n",
    "```\n",
    "where data is the attribute name.\n",
    "\n",
    "You acceess object methods with the syntax\n",
    "```python\n",
    "object.method(arguments)\n",
    "```\n",
    "\n",
    "For example, for a list \"`append` is a method that was covered in the introductory lesson.\n",
    "\n",
    "```\n",
    "my_list = []\n",
    "my_list.append(1) # \"append\" is a method\n",
    "```\n",
    "</div>    \n",
    "\n",
    "In this lesson, we will create and manipulate RDKit `mol` objects. RDKit `mol` objects represent molecules and have\n",
    "attributes (data) and methods (actions) associated with molecules.\n",
    "\n",
    "We are going to use a part of RDKit called `Chem`. To use `Chem`, we first have to import it. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afc9bc63",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from rdkit import Chem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5b835ad",
   "metadata": {},
   "source": [
    "## Creating Molecules with RDKit\n",
    "\n",
    "To get information about molecules in RDKit, we have to first create objects representing molecules. RDKit has molecule object that can be used to retrieve information or calculate properties. To create a molecule object, we have to communicate the molecule identity in a way that computers understand. \n",
    "\n",
    "### Creating molecules using SMILES\n",
    "\n",
    "In the last lesson, we learned about molecular representations using SMILES strings. Now we will use SMILES strings to create molecule objects in RDKIT. \n",
    "\n",
    "We can create a representation of methane using RDKit by using the `MolFromSmiles` function in `rdkit.Chem`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40b29a7c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "methane = Chem.MolFromSmiles(\"C\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da5b102b-e997-47b4-8a0d-7cc5befc1b7e",
   "metadata": {},
   "source": [
    "The `methane` variable is now an RDKit \"molecule object\".\n",
    "\n",
    "In a Jupyter environment, putting a variable that is an RDKit molecule as the only or last thing in a cell,\n",
    "will result in a picture of the molecule as an output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99734c44",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "methane"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e98c39d3",
   "metadata": {},
   "source": [
    "<div class=\"exercise admonition\">\n",
    "<p class=\"admonition-title\">Check Your Understanding</p>\n",
    "<p> Create RDKit molecules for the following molecules. You can look up the SMILES strings on <a href=\"https://pubchem.ncbi.nlm.nih.gov/\">PubChem</a>:\n",
    "<p>\n",
    "    <ul>\n",
    "        <li> Propane\n",
    "        <li> Ethene\n",
    "        <li> Cyclohexane\n",
    "        <li> Benzene\n",
    "    </ul>\n",
    "</p>\n",
    "<p>Create variables for each molecule. The variable names should be the molecule name (all lowercase)</p>\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bbae7a92-6647-4361-b5ed-f45363f4d249",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Your answers go here\n",
    "propane = \n",
    "ethene = \n",
    "cyclohexane = \n",
    "benzene = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "395bc20c-bef9-4941-a2f8-fca830baae2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try visualizing your molecules here - you visualize a molecule as putting its variable name as the last thing in a cell\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e22a13a2-305e-458b-9caa-ba636b836e1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try visualizing your molecules here - you visualize a molecule as putting its variable name as the last thing in a cell\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fc9b578-6889-43d9-85c0-0f07111053b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try visualizing your molecules here - you visualize a molecule as putting its variable name as the last thing in a cell\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0942cfe3-0ca2-4145-baaf-c7dff83fed24",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try visualizing your molecules here - you visualize a molecule as putting its variable name as the last thing in a cell\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95ad0cdb",
   "metadata": {},
   "source": [
    "## RDKit and Molecular Properties\n",
    "\n",
    "RDKit will allow us to access molecular properties of the molecules we've created.\n",
    "For example, we can have RDKit tell us the number of atoms in a molecule using the `GetNumAtoms` methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e045598b-847c-4d21-a336-f4b6cd7bf0b0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "num_methane = methane.GetNumAtoms()\n",
    "print(f\"The number of atoms in methane is {num_methane}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc5d6ee6",
   "metadata": {},
   "source": [
    "This isn't exactly what we expect. By default, RDKit only counts \"heavy atoms\". This means that hydrogen isn't included. We can tell RDKit to count hydrogens by adding `onlyExplicit=False` to our `GetNumAtoms` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e94e9f4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "num_methane_h = methane.GetNumAtoms(onlyExplicit=False)\n",
    "print(f\"The number of atoms in methane including hydrogens is {num_methane_h}.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e78625d",
   "metadata": {},
   "source": [
    "We can also get information about atoms or  bonds from RDKit molecules. We can use the `GetAtoms()` function, then use a loop to look at t`for` loop to look at the atoms. In the `for` loop, we can get information about each atom, like the atomic symbol or mass."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bf913ea",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ethanol = Chem.MolFromSmiles(\"CCO\")\n",
    "\n",
    "for atom in ethanol.GetAtoms():\n",
    "    print(atom.GetSymbol(), atom.GetMass())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e14d40c2",
   "metadata": {},
   "source": [
    "Like the previous example, the molecule does not have hydrogens by default.\n",
    "If you need to work with the hydrogens, you can use a function to add them called `AddHs`.\n",
    "Now, when we look through the atoms, we get the hydrogens as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "107fe76e",
   "metadata": {},
   "outputs": [],
   "source": [
    "ethanol_h = Chem.AddHs(ethanol)\n",
    "\n",
    "ethanol_h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a96a43e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "for atom in ethanol_h.GetAtoms():\n",
    "    print(atom.GetSymbol())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6dc6fb9",
   "metadata": {},
   "source": [
    "You can use the `GetBonds` method to get the bonds in a molecule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "371244ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "for bond in ethanol.GetBonds():\n",
    "    print(bond.GetBondType())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e9469fa-d7c2-44e6-b83c-6f06946503ba",
   "metadata": {},
   "source": [
    "<div class=\"exercise admonition\">\n",
    "<p class=\"admonition-title\">Check Your Understanding</p>\n",
    "<p> Print atom and bond info for one or more of the molecules you created in the previous exercise. </p>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7027bb2e-b388-4c5f-a055-27d989641d25",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec52c9cd-a5fa-44f6-9e68-f5b2086df591",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "b39599a7",
   "metadata": {},
   "source": [
    "## RDKit and Molecular Descriptors\n",
    "\n",
    "A molecular descriptor is a numerical value or a set of values that represent specific structural or chemical features of a molecule. Molecular descriptors are based on molecular structure and allow statistical analysis and similarity measurements on molecules.\n",
    "\n",
    "RDKit supports the calculation of many molecular descriptors. You can see a [full list of RDKit descriptors](https://www.rdkit.org/docs/GettingStartedInPython.html#list-of-available-descriptors) or [see the module documentation](https://www.rdkit.org/docs/source/rdkit.Chem.Descriptors.html).\n",
    "\n",
    "To get molecular descriptors from RDKit, we import the `Descriptors` module.\n",
    "\n",
    "```python\n",
    "from rdkit.Chem import Descriptors\n",
    "```\n",
    "\n",
    "To get a descriptor, you do\n",
    "\n",
    "```python\n",
    "Descriptors.descriptor_name(molecule_variable)\n",
    "\n",
    "```\n",
    "\n",
    "RDKit supports a number of molecular descriptors, a list of which is given below, along with the RDKit method for the property.\n",
    "\n",
    "\n",
    "Name of Property      | Name of RDKit Descriptor Function\n",
    "----------------------|-----------------------------------|\n",
    "molecular weight      | Descriptors.MolWt\n",
    "number of heavy atoms | Descriptors.HeavyAtomCount\n",
    "number of H-bond donors| Descriptors.NumHDonors\n",
    "number of H-bond donors| Descriptors.NumHAcceptors\n",
    "octanol-water partition coefficient| Descriptors.MolLogP\n",
    "topological polar surface area | Descriptors.TPSA\n",
    "number of rotatable bonds      | Descriptors.NumRotatableBonds\n",
    "number of aromatic rings       | Descriptors.NumAromaticRings\n",
    "number of aliphatic rings      | Descriptors.NumAliphaticRings\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d8d3ac3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from rdkit.Chem import Descriptors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3ed7282",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(\"Printing info for methane:\")\n",
    "print(f\"The molecular weight is {Descriptors.MolWt(methane)}\")\n",
    "print(f\"The number of aromatic rings is {Descriptors.NumAromaticRings(methane)}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02f13392-9307-4398-8f31-c23b6f5faa60",
   "metadata": {},
   "source": [
    "<div class=\"exercise admonition\">\n",
    "<p class=\"admonition-title\">Check Your Understanding</p>\n",
    "    <p> Print molecular descriptors for some of the molecules you've created</p>\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb3aa2de-f22c-4408-addb-11a1417022eb",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "401eb6c2-179f-4e91-95cf-b9799dc66108",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6868227-ea3c-48dc-8759-1cd8db8830e2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5b9f1e8-469c-42db-b628-03b4c0bb0c1f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4143ac73-a5b6-419d-aa8f-692b33bb83ea",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "8a5be400",
   "metadata": {},
   "source": [
    "## RDKit and Stereochemistry\n",
    "\n",
    "RDKit can be used to find chiral centers and print stereochemistry information about molecules. Consider 2-chlorobutane, created below.\n",
    "We are using a canonical SMILES which does not include stereo center information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f05c1617-9d2f-4f83-b151-2d15aa8168c3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "chlorobutane = Chem.MolFromSmiles(\"CCC(C)Cl\")\n",
    "chlorobutane"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60375660-0e25-49ef-aa16-1c89637018f9",
   "metadata": {},
   "source": [
    "However, RDKit can analyze the molecule and tell us about the number of chiral centers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b091f596",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Can tell us how many stereocenters\n",
    "Chem.rdMolDescriptors.CalcNumAtomStereoCenters(chlorobutane)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d9c32b8-0ab1-4348-a5c3-6d81144706f7",
   "metadata": {},
   "source": [
    "We can also retrieve the atom indices of the chiral centers and see their orientation.\n",
    "If the orientation was not specified, RDKit will list it as `?`. \n",
    "We also have to tell RDKit to tell us about unassigned chiral centers using\n",
    "`includeUnassigned=True`.\n",
    "\"Unassigned\" refers to the fact that we have not specified the stereochemistry.\n",
    "\n",
    "We can get the atom index of the stereocenter using `FindMolChiralCenters`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e2c4e5a-4f1b-4e88-aeaf-6b3585088e2b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Can tell us the atom number and the orientation.\n",
    "Chem.FindMolChiralCenters(chlorobutane, includeUnassigned=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea1457c4-0a0c-439f-aa99-a07e8cfea28e",
   "metadata": {},
   "source": [
    "### Visualizing Chiral Centers\n",
    "\n",
    "Although we can visualize molecules by making them the last thing in a cell, RDKit also allows us to make more complicated visualizations using its\n",
    "`Draw` module. To make the same visualization we did earlier we could do the following"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "decd78e2-c887-4240-a0cd-d77e965641fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit.Chem import Draw\n",
    "\n",
    "Draw.MolToImage(chlorobutane)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff5203a6-338b-470c-82aa-3631551a4614",
   "metadata": {},
   "source": [
    "The `MolToImage` function can be used to highlight particular atoms. \n",
    "To highlight the chiral center in this molecule, we might first use\n",
    "`Chem.FindMolChiralCenters` to get the index of the chiral centers,\n",
    "then pass the atom indices of the chiral centers into `MolToImage`, as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd9b2525-7250-4b11-a5bf-f862968f4b7a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "chiral_centers = Chem.FindMolChiralCenters(chlorobutane, includeUnassigned=True)\n",
    "print(chiral_centers)\n",
    "\n",
    "# We use 2 because that is the index of the chiral center.\n",
    "Draw.MolToImage(chlorobutane, highlightAtoms=[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da6a8dff-6d9d-4439-aca7-984e98fc365f",
   "metadata": {},
   "source": [
    "We could have also used Python instead of filling in the number.\n",
    "A `for` loop would make this more transferable between different atoms (hint for the next exercise!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a19d0e2f-1d38-40e2-a3b2-cb69be0480bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "chiral_centers = Chem.FindMolChiralCenters(chlorobutane, includeUnassigned=True)\n",
    "\n",
    "chiral_indices = []\n",
    "\n",
    "for center in chiral_centers:\n",
    "    chiral_indices.append(center[0])\n",
    "\n",
    "Draw.MolToImage(chlorobutane, highlightAtoms=chiral_indices)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9c43cdd-b9ca-49e1-8fd3-520bb4ed4aa2",
   "metadata": {},
   "source": [
    "<div class=\"exercise admonition\">\n",
    "<p class=\"admonition-title\">Exercise</p>\n",
    "    <p> Create a molecule using the SMILES string C1[C@@H](Cl)[C@H](Cl)CCC1. Print the stereochemistry information and\n",
    "        create a visualization highlighting the stereocenters.</p>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23991157-3516-448f-8159-9d777beab358",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a398d681-9a75-4b7f-99d2-9dc60b74730b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f3019c1-fb23-40d5-953d-619534e17448",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31444b0f-a345-4b58-b720-6fcd6166f6bb",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "e976bbd5-a457-46cc-8196-d6da501b2aaf",
   "metadata": {},
   "source": [
    "## RDKit and Molecular Similarity\n",
    "\n",
    "### Substructure Searches\n",
    "\n",
    "RDKit can be used to search molecules for substructures. \n",
    "You can define your own substructures using SMILES. \n",
    "For example, we'll first look for carbons in caffeine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6da72486-c3b4-4bbe-9613-2985090d8b91",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "caffeine = Chem.MolFromSmiles('CN1C=NC2=C1C(=O)N(C(=O)N2C)C')\n",
    "carbon = Chem.MolFromSmiles('C')\n",
    "\n",
    "matches = caffeine.GetSubstructMatches(carbon)\n",
    "\n",
    "caffeine"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3027482-0ffd-4124-9a07-d7a9fde69f77",
   "metadata": {},
   "source": [
    "In this example, you can see that all of the carbons are highlighted. \n",
    "\n",
    "The code `matches = caffeine.GetSubstructMatches(carbon)` will highlight all of the carbons in the caffeine molecule because it is searching for a specific substructure, which in this case is just a single carbon atom.\n",
    "\n",
    "But, what if we just wanted to highlight the methyl groups? We would have to use a different strategy for substructure matching.\n",
    "\n",
    "#### SMARTS Strings\n",
    "A [SMARTS (SMiles ARbitrary Target Specification)](https://en.wikipedia.org/wiki/SMILES_arbitrary_target_specification) string is a pattern-matching language for specifying substructures in molecules. It is similar to SMILES but has additional features for specifying atom and bond properties. In contrast to SMILES, SMARTS allows for specifying substructures based on patterns of atoms and bonds, rather than just specifying the exact arrangement of atoms and bonds in a molecule.\n",
    "\n",
    "In this case, we want to highlight the methyl group in caffeine, which contains a single carbon atom bonded to three hydrogen atoms. We can use a SMARTS string to define this substructure as `[CH3]`, which represents a carbon atom with three hydrogen atoms bonded to it.\n",
    "\n",
    "Using a SMARTS string is necessary in this case because simply searching for a single carbon atom with the code `carbon = Chem.MolFromSmiles('C')` would highlight all carbons in the molecule, not just the methyl group. By using a SMARTS string, we can specify the exact substructure we want to highlight and avoid highlighting unintended parts of the molecule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65d9212b-1b66-4c23-8314-2b9b10c7b817",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "caffeine = Chem.MolFromSmiles('CN1C=NC2=C1C(=O)N(C(=O)N2C)C')\n",
    "methyl_pattern = Chem.MolFromSmarts('[CH3]')\n",
    "\n",
    "# Get the indices of the matching atoms\n",
    "matches = caffeine.GetSubstructMatches(methyl_pattern)\n",
    "caffeine"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4d29c59-6587-4ad1-adba-907ff11c631a",
   "metadata": {},
   "source": [
    "You can see a list of SMARTS strings for different functional groups [at this link](https://www.daylight.com/dayhtml_tutorials/languages/smarts/smarts_examples.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4acb103b-f24b-4276-901b-a6f2e35e180b",
   "metadata": {},
   "source": [
    "<div class=\"exercise admonition\">\n",
    "<p class=\"admonition-title\">Exercise</p>\n",
    "    <p>Create an RDKit molecule object for aspirin (acetylsalicylic acid) using its SMILES string \"CC(=O)OC1=CC=CC=C1C(=O)O\".</p>\n",
    "    <p>Visualize the molecule, highlighting any carboxyl groups (COOH) in the structure. Use a SMARTS string to define the carboxyl group substructure.</p>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cde5bf62-668a-4141-a6b1-afb7a44e8ed8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8aa2579f-2598-47a8-a277-c4274bd03c30",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "dc1b8369-d480-4a64-8b57-04eedece3244",
   "metadata": {},
   "source": [
    "## Molecular Fingerprints\n",
    "\n",
    "Molecular fingerprinting is a method that represents the structure of a molecule in a way so that molecules may be compared. There are numerous molecular fingerprint algorithms that you might choose, each with their own strengths and weaknesses. You can see more information on fingerprinting in RDKit in [this presentation](https://www.rdkit.org/UGM/2012/Landrum_RDKit_UGM.Fingerprints.Final.pptx.pdf) or in the [RDKit documentation](https://www.rdkit.org/docs/GettingStartedInPython.html#fingerprinting-and-molecular-similarity).\n",
    "\n",
    "Let's see how we might use RDKit to compare some molecules using fingerprints.\n",
    "We already have aspirin defined. We will also create a molecule for salicylic acid and one for acetaminophen.\n",
    "Aspirin is derived from salicylic acid, so we would expect the structures to be similar.\n",
    "Acetaminophen is another molecule that is used as a pain reliever."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "124b129d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit import Chem\n",
    "from rdkit import DataStructs\n",
    "from rdkit.Chem import Draw\n",
    "\n",
    "salicylic_acid_smiles = \" C1=CC=C(C(=C1)C(=O)O)O\"\n",
    "aspirin_smiles = \"CC(=O)OC1=CC=CC=C1C(=O)O\"\n",
    "acetaminophen_smiles = \"CC(=O)NC1=CC=C(C=C1)O\"\n",
    "\n",
    "aspirin = Chem.MolFromSmiles(aspirin_smiles)\n",
    "salicylic_acid = Chem.MolFromSmiles(salicylic_acid_smiles)\n",
    "acetaminophen = Chem.MolFromSmiles(acetaminophen_smiles)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c35dc81-ea82-471c-9d32-2b1480ed4dbc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "salicylic_acid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3b49b61-4b82-4c52-9b6f-a87b96dd48ea",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "aspirin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9c7e8c9-9069-4097-b09e-f7efbd6f2f21",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "acetaminophen"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e56c6c8-1163-4005-8b4a-d5c48b6509bf",
   "metadata": {},
   "source": [
    "From inspecting these molecules visually, we can see that aspirin and salicylic acid look more similar than aspirin and acetaminophen. \n",
    "We will see this next using a molecular fingerprint.\n",
    "\n",
    "We will use the `Chem.RDKFingerprint`, which is an RDKit specific fingerprint, though there are other options."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "864cf1eb-7eca-4b71-a490-53078492a83c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "acetaminophen_fp = Chem.RDKFingerprint(acetaminophen)\n",
    "aspirin_fp = Chem.RDKFingerprint(aspirin)\n",
    "salicylic_acid_fp = Chem.RDKFingerprint(salicylic_acid)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90af75ae-ae6c-4513-842f-3e3ec76cb033",
   "metadata": {},
   "source": [
    "To calculate the similarity between molecular fingerprints, a measure called the Tanimoto similarity is often used. \n",
    "This number ranges from 0 to 1.\n",
    "A higher number means that molecules are more similar. For example, the similarity of a molecule with itself would be 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "927f103e-b65e-464c-b5b4-54b641b92fa2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "similarity = DataStructs.TanimotoSimilarity(aspirin_fp, salicylic_acid_fp)\n",
    "print(similarity)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2792e17-a2b3-430c-8dd1-b79d283bd5b4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "similarity = DataStructs.TanimotoSimilarity(aspirin_fp, acetaminophen_fp)\n",
    "print(similarity)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc73d27a-f9ea-4595-a5db-66704918354e",
   "metadata": {},
   "source": [
    "As expected, our similarity scores indicate that aspirin and salicylic acid are more similar than aspirin and acetaminophen. \n",
    "This starts being more useful when working with larger data sets."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28a48040-bb26-44b1-b119-c51b9de140b4",
   "metadata": {
    "tags": []
   },
   "source": [
    "<div class=\"exercise admonition\">\n",
    "<p class=\"admonition-title\">Final Challenge</p>\n",
    "\n",
    "<p> For the following sets of molecules, generate the RDKit molecule using a SMILES string (you can find this on PubChem).\n",
    "    Identify which molecules contain an amine group and which contain a phenyl group. Finally, use fingerprints to compute\n",
    "    the similarity between a few of the molecules.</p>\n",
    "\n",
    "<ul>\n",
    "<li>Neurotransmitters:</li>\n",
    "    <ul>\n",
    "        <li> Dopamine </li>\n",
    "        <li> Serotonin </li>\n",
    "        <li> Norephinephrine </li>\n",
    "    </ul>\n",
    "<li>Hormones:</li>\n",
    "    <ul>\n",
    "        <li> Cortisol </li>\n",
    "        <li> Testosterone </li>\n",
    "        <li> Estradiol </li>\n",
    "    </ul>\n",
    "<ul>\n",
    "\n",
    "</div>\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c10a9ee8-f382-4f2b-a98f-4aadf16959b2",
   "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.10.12"
  },
  "vscode": {
   "interpreter": {
    "hash": "16d4a7bb199d969b1271ebe46f77414b0d9cd01b3c3983c2b2742fc6cd4503d3"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}