03_molecular_similarity.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "10110027-0410-40ea-8950-8e4047699603",
   "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",
    "# Cheminformatics: Molecular Similarity and Molecular Descriptors\n",
    "\n",
    "<strong>Author(s):</strong> Jessica A. Nash, The Molecular Sciences Software Institute\n",
    "\n",
    "<div class=\"alert alert-block alert-info\"> \n",
    "<h2>Overview</h2>\n",
    "\n",
    "<strong>Questions:</strong>\n",
    "\n",
    "* How are molecules represented on computers for cheminformatics applications?\n",
    "\n",
    "* How can I measure the similarity of two molecules?\n",
    "\n",
    "* How can I find out if a molecule contains a particular substructure?\n",
    "\n",
    "<strong>Objectives:</strong>\n",
    "\n",
    "* Learn about graph representation of molecules.\n",
    "\n",
    "* Learn about molecular fingerprints.\n",
    "\n",
    "* Use RDKit to measure molecular similarity.\n",
    "\n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">\n",
    "<h3>Need the solutions?</h3>\n",
    "\n",
    "This notebok contains \"fill-in-the-blank\" cells and other exercises. If you would like to access the solutions, please see\n",
    "[this notebook](filled_notebooks/03_molecular_similarity.ipynb).\n",
    "\n",
    "</div>\n",
    "\n",
    "\n",
    "Cheminformatics is often used in the context of molecular discovery. \n",
    "For example, you might want to process a large number of molecules to find molecules that are similar to a target, or that contain a particular functional group. \n",
    "But how can you tell if two molecules are similar?\n",
    "\n",
    "A chemist looking at two molecules might be able to make a judgment based on their chemical knowledge. However, when dealing with a large number of molecules, this problem can become intractable. There are ways to detect molecular similarity on a computer that allow fast comparison and processing of molecules. \n",
    "\n",
    "## Graph Representation of Molecules\n",
    "Graph theory is a branch of mathematics that studies the properties and applications of graphs, which are structures made up of nodes (or vertices) connected by edges (or lines).\n",
    "A molecule can be represented as a graph where atoms are nodes and bonds are edges. This representation allows the use of various graph-based algorithms to analyze molecular structures. \n",
    "The figure below (borrowed from Wikipedia) shows a graph with numbered nodes (circles) connected by edges (lines).\n",
    "\n",
    "\n",
    "<div style=\"text-align:center;\">\n",
    "  <img src=\"images/molecular_graph.png\" style=\"display: block; margin: 0 auto; max-height:300px;\">\n",
    "</div>\n",
    "\n",
    "Image from [Wikipedia: Graph theory](https://en.wikipedia.org/wiki/Graph_theory)\n",
    "\n",
    "As a chemist, you probably don't have a hard time imagining a molecule as a graph. The atoms would be \"nodes\" in the graph, while the bonds would be \"edges\".\n",
    "One can also represent bond order by changing the \"weight\" of graph edges.\n",
    "In this view, one way to tell if molecules are the same is to check that their molecular graphs are the same.\n",
    "\n",
    "### Substructure Searches\n",
    "\n",
    "A substructure search is a cheminformatics technique used to identify molecules that contain a specific pattern or structure within a larger molecule.\n",
    "\n",
    "The figure below shows the beneze substructure matched in the aspirin molecule.\n",
    "\n",
    "<div style=\"text-align:center;\">\n",
    "  <img src=\"images/substructure.png\" style=\"display: block; margin: 0 auto; max-height:300px;\">\n",
    "</div>\n",
    "\n",
    "\n",
    "One way this can be done is with graph representaiton and graph theory. \n",
    "When doing substructure search using RDKit, it can check to see if the molecular graph for one molecule contains the molecular graph of a smaller molecule.\n",
    "This is a concept called \"subgraph isomorphism\". You can see more about this topic in this [presentation](https://www.rdkit.org/UGM/2012/Landrum_RDKit_UGM.Fingerprints.Final.pptx.pdf).\n",
    "\n",
    "This is all very nice information, but you actually don't have to fully understand it in order to do a substructure search.\n",
    "In the cells below, we demonstrate searching molecules for a smaller structure using substructure search. \n",
    "As you can see from the cells, this functionality is very nicely built into RDKit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a58f4581-13c4-478d-bafe-1e0846d3b412",
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit import Chem\n",
    "\n",
    "caffeine = Chem.MolFromSmiles('CN1C=NC2=C1C(=O)N(C(=O)N2C)C')\n",
    "carbonyl = Chem.MolFromSmiles('C=O')\n",
    "\n",
    "matches = caffeine.GetSubstructMatches(carbonyl)\n",
    "\n",
    "caffeine"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32a4c55b-362f-472b-a784-105b5c4fe015",
   "metadata": {},
   "source": [
    "As you can see above, we can define two molecules and search for one inside of the other. \n",
    "\n",
    "<div class=\"alert alert-block alert-warning\"> \n",
    "<h3>Check Your Understanding</h3>\n",
    "\n",
    "Using the skills you've learned already, perform a substructure search for methyl groups in caffeine in the cell below.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efabcc57-46cf-4742-b39a-e0261307c07f",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd03244e-945a-4265-b544-10a2dbffd5e4",
   "metadata": {},
   "source": [
    "From the substructure search in the cell above, you are probably seeing that you're not matching exactly what you want. In order to find just the methyl groups, we have to use yet another molecular representation.\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')` highlights 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": "1b5bba4e-91de-4c91-986b-049e0cac168c",
   "metadata": {},
   "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": "89d4a8b5-0454-46bd-9805-8e72a60f08ef",
   "metadata": {},
   "source": [
    "\n",
    "<div class=\"alert alert-block alert-warning\"> \n",
    "<h3>Exercise</h3>\n",
    "\n",
    "<p>Create an RDKit molecule object for aspirin. You can find its SMILES representation by using PubChem.</p>\n",
    "<p>Visualize the molecule, highlighting any carboxyl groups in the structure. Use a SMARTS string to define the carboxyl group substructure.</p>\n",
    "\n",
    "<strong>Important:</strong> Unless you know the rules for writing SMARTS strings (which you would have learned somewhere else, because we didn't cover it), it's unlikely that you'll guess what the SMARTS for carboxyl is. 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).\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03bf284b-6730-43b6-8c5b-242d095853eb",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "4a25724d-2e54-4a1a-b5b8-3a0da72b8e97",
   "metadata": {},
   "source": [
    "## Molecular Fingerprints\n",
    "Molecular fingerprints are representations of molecules that are usually bit strings, or vectors of 0's and 1's. Fingerprints are built by considering the molecular structure (often as a graph representation) and applying a certain algorithm to create the vector. The image below shows a representation of a molecular fingerprint. The molecule is separated into different parts with each part setting a bit (changing a 0 to a 1) in the fingerprint.\n",
    "\n",
    "![molecular fingerprint](images/molecular_fingerprint.png)\n",
    "\n",
    "Image from [Chemistry LibreTexts](https://chem.libretexts.org/Courses/Intercollegiate_Courses/Cheminformatics/06%3A_Molecular_Similarity/6.01%3A_Molecular_Descriptors) [Cheminformatics Course](https://chem.libretexts.org/Courses/Intercollegiate_Courses/Cheminformatics)\n",
    "\n",
    "The exact fingerprint will depend on the algorithm that is used to construct it.\n",
    "There are many different fingerprinting algorithms. But they tend to fall into two types of groups - similarity or substructure fingerprints. A common similarity fingerprint that is used is the Morgan fingerprint. A common substructure fingerprint that is used is the [Daylight fingerprint](https://www.daylight.com/dayhtml/doc/theory/theory.finger.html) (the RDKFingerprint is a Daylight-like fingerprint).\n",
    "\n",
    "The cell below calculates the RDKit fingerprint for caffeine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca97af7e-cc2c-49e1-ab71-7e8058106131",
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit.Chem import AllChem\n",
    "\n",
    "aspirin = Chem.MolFromSmiles(\"CC(=O)Oc1ccccc1C(=O)O\")\n",
    "benzene = Chem.MolFromSmiles(\"c1ccccc1\")\n",
    "\n",
    "aspirin_fingerprint = AllChem.RDKFingerprint(aspirin)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3ef1152-d5c1-40a9-a66d-0a0f3c83facc",
   "metadata": {},
   "source": [
    "You usually wouldn't interact with a fingerprint directly, but using RDKit we can examine it to understand more about how it works.\n",
    "We can print the bit string representation using `.ToBitString`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71d91ca4-4848-46ee-9441-9bb81bec48ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "aspirin_fingerprint.ToBitString()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c79828c2-5b84-43f0-bb2d-f04d95141d9a",
   "metadata": {},
   "source": [
    "If you print the benzene string, you will see a string of the same length, but with fewer bits set.\n",
    "There are fewer bits because there are fewer molecular patterns in the smaller molecule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9018a8f0-0e32-4e32-8cad-377e5cfe5838",
   "metadata": {},
   "outputs": [],
   "source": [
    "benzene_fingerprint = AllChem.RDKFingerprint(benzene)\n",
    "benzene_string = benzene_fingerprint.ToBitString()\n",
    "\n",
    "print(benzene_string)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "781dd591-b789-4283-93b2-ffe684ee61c9",
   "metadata": {},
   "source": [
    "Using this bit string, we can do a number of different analyses. \n",
    "For example, substructure searches using bit strings check that bits that are set in the smaller molecule are also set in the larger target molecule.\n",
    "\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\"> \n",
    "<h3>Exercise</h3>\n",
    "\n",
    "Fill in the for loop below to check that bits set in the aspiring fingerprint string are also set in the benzene fingerprint string.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0869b43-18dd-4714-af92-1de49932eb2b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "9778634e-0344-4d87-a406-c0ab48b8ff75",
   "metadata": {},
   "source": [
    "The above output tells us that every bit that is set in the benzene fingerprint is also set in the aspirin target. Fingerprint substructure searches are often used for faster, but less accurate substructure searches, and are commonly used for screen large numbers of molecules quickly."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9e5ecb4-5348-4f3d-b519-4fbff21f666a",
   "metadata": {},
   "source": [
    "### Measuring Similarity\n",
    "\n",
    "Because the fingerprints are numbers, we can measure how similar two fingerprints are using different similarity metrics.\n",
    "One common similarity metric is called the Tanimoto similarity.\n",
    " The Tanimoto similarity is calculated as follows:\n",
    "\n",
    "$$\n",
    "T(A, B) = \\frac{A \\cap B}{A + B - A \\cap B}\n",
    "$$\n",
    "\n",
    "Where:\n",
    "- $A$ and $B$: are the sets of bits in the fingerprint vectors for molecules $A$ and $B$ respectively.\n",
    "- $A \\cap B$: This represents the intersection of sets $A$ and $B$, i.e., the number of bits that are '1' (set) in both $A$ and $B$.\n",
    "- $A + B$: This is the sum of all '1' bits in both $A$ and $B$.\n",
    "- $A + B - A \\cap B$: This term represents the union of sets $A$ and $B$, calculated as the total number of unique '1' bits across both fingerprints.\n",
    "\n",
    "The Tanimoto similarity ranges from 0.0 to 1.0, with 1.0 representing identical fingerprints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "820dd583-c661-4db5-b510-aa204e0c0f26",
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit import DataStructs\n",
    "\n",
    "DataStructs.TanimotoSimilarity(benzene_fingerprint, aspirin_fingerprint)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75edf094-0dbb-4094-93bc-cf0cc113c51d",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\"> \n",
    "<h3>Exercise</h3>\n",
    "\n",
    "Use the Tanimoto similarity to measure if benzene is more similar to pyridine or aniline using the RDKit fingerprint.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a47675e2-4294-4625-a71a-4b411566f908",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "92228613-2c40-46c0-823b-63ec2700da07",
   "metadata": {},
   "source": [
    "## Molecular Descriptors\n",
    "\n",
    "The fingerprints discussed above are examples of **molecular descriptors**.\n",
    "\n",
    "A molecular descriptor is a numerical value that represents some property of a molecule.\n",
    "If we want to be able to create a model that describes chemical behavior, we have to be able to convert information about molecules into numerical representations.\n",
    "This is where the concept of a molecular descriptor comes in. \n",
    "\n",
    "Descriptors can be 0 dimensional (molecular weight, number of heavy atoms, etc.), 1 dimensional (counts of atom types, hydrogen bond donors/acceptors), 2 dimensional (fingerprints, other graph representations), 3 dimensional (polar surface area). The dimensionality of the descriptor defines what kind of dimensional information you need about the molecule in order to define the descriptor. For example, the fingerprint, a 2D descriptor, depends on the connectivity, or 2 dimensional structure of the molecule.\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",
    "For example, we can calculate the molecular weight of our caffeine molecule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e860c2f-672b-419f-ab2c-ee58602ab422",
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit.Chem import Descriptors\n",
    "\n",
    "caffeine_mol_wt = Descriptors.MolWt(caffeine)\n",
    "print(caffeine_mol_wt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7986dcba-1058-49f6-9788-5433288377bd",
   "metadata": {},
   "source": [
    "RDKit provides hundreds of molecular descriptors.\n",
    "You can calculate all of the molecular descriptors for a molecule at once using `CalcMolDescriptors`.\n",
    "\n",
    "As an example of how we can apply molecular descriptors to make predictions, we will consider Lipinski's Rule of 5.\n",
    "Lipinski's Rule of 5 is a **guideline** that helps determine if a drug is likely to be absorbed well by the body. It states that good oral drugs typically have no more than 5 hydrogen bond donors, 10 hydrogen bond acceptors, a molecular weight under 500 daltons, and a log P (measure of solubility) under 5.\n",
    "\n",
    "* Molecular Weight <= 500 Da\n",
    "* No. Hydrogen Bond Donors <= 5\n",
    "* No. Hydrogen Bond Acceptors <= 10\n",
    "* LogP <= 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "988dcfa9-389a-42c4-bb9b-4f68666f4583",
   "metadata": {},
   "outputs": [],
   "source": [
    "MW = Descriptors.MolWt(aspirin)\n",
    "HBA = Descriptors.NOCount(aspirin)\n",
    "HBD = Descriptors.NHOHCount(aspirin)\n",
    "LogP = Descriptors.MolLogP(aspirin)\n",
    "\n",
    "rules = [ MW <= 500, HBD <=5, HBA <=10, LogP <=5 ]\n",
    "print(rules)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cc4c878-c5c6-4cbc-9e10-91b9760b9201",
   "metadata": {},
   "source": [
    "Based on these rules, aspirin is expected to be well-absorbed by the body.\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\"> \n",
    "<h2>Final Challenge</h2>\n",
    "\n",
    "Check Lipinski's Rule of 5 for [insulin](https://pubchem.ncbi.nlm.nih.gov/compound/Insulin#section=Canonical-SMILES).\n",
    "Do you expect that insulin will be well-absorbed when taken orally based on these criteria?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c0aa9ee-1f21-40c6-ae04-86c6abdeb245",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f00ab72a-2e36-4562-b5ad-e829a2a93e8b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c22fe39-8bb3-4476-a11f-56e4e8f6508d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python3.10 (MolSSI)",
   "language": "python",
   "name": "molssi"
  },
  "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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}