{ "cells": [ { "cell_type": "code", "execution_count": null, "id": "infrared-polish", "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "from rdkit import Chem\n", "from rdkit.Chem import AllChem\n", "import nglview\n", "import csv\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "id": "7042051a", "metadata": { "tags": [] }, "source": [ "### As a starting point, you have been given a list of molecules in SMILES format containing compounds taken from the WHO (World Health Organization) Essential Medicine List. \n", "\n", "Many, but not all of the essential medicines are small molecules. Roughly speaking, a small molecule is a molecule of relatively low molecular weight - a commonly cited value is < 900 daltons - in order to distinguish them from larger molecules such as proteins and polymers.\n", "\n", "`WHO_EML_SMILES.csv` is a text file, where on each line there are two values separated by a comma. The first value is the name of the drug and the second value contains the 2D molecular structure in SMILES format. These SMILES strings were generated as part of a [research paper](https://doi.org/10.1039/C9RE00348G) and kindly provided by the authors, Klavs Jensen and Hanyu Gao.\n", "\n", "We will begin by using the `csv` module to read in the file. As we can see here, the csv.reader function returns an iterable object containing the contents of the file." ] }, { "cell_type": "code", "execution_count": null, "id": "376ee787", "metadata": {}, "outputs": [], "source": [ "with open('WHO_EML_SMILES.csv') as csv_file:\n", " reader = csv.reader(csv_file)\n", " for line in reader:\n", " print(line)" ] }, { "cell_type": "markdown", "id": "9912b4fb", "metadata": {}, "source": [ "Take abacavir for example, the first line in the file. This is a medication used to treat HIV/AIDS. To start working with this molecule, we could build the structure in Avogadro, but a more systematic way would be to create a molecule object directly from the SMILES string. The RDKit toolkit allows us to do this." ] }, { "cell_type": "code", "execution_count": null, "id": "static-expert", "metadata": {}, "outputs": [], "source": [ "# Generate an RDKit molecule object from a SMILES string.\n", "# The example uses abacavir taken from the above printout.\n", "m1 = Chem.MolFromSmiles('Nc1nc(NC2CC2)c3ncn([C@@H]4C[C@H](CO)C=C4)c3n1')" ] }, { "cell_type": "code", "execution_count": null, "id": "parallel-symbol", "metadata": {}, "outputs": [], "source": [ "# This is how the RDKit molecule object is represented - as a 2D drawing.\n", "m1" ] }, { "cell_type": "code", "execution_count": null, "id": "opposed-presentation", "metadata": {}, "outputs": [], "source": [ "# Starting from SMILES, there are no 3-D coordinates and they need to be generated.\n", "# This is how RDKit generates 3-D conformers for a molecule.\n", "# Under the hood, it uses a method called ETDKG to do this.\n", "AllChem.EmbedMolecule(m1)" ] }, { "cell_type": "code", "execution_count": null, "id": "looking-nebraska", "metadata": {}, "outputs": [], "source": [ "# Now we can print out the 3D coordinates. Note that hydrogen atoms are missing.\n", "print(Chem.MolToXYZBlock(m1))" ] }, { "cell_type": "code", "execution_count": null, "id": "settled-coral", "metadata": {}, "outputs": [], "source": [ "# nglview has a nice tool for allowing us to view the 3D coordinates as well.\n", "nglview.show_rdkit(m1)" ] }, { "cell_type": "code", "execution_count": null, "id": "thousand-diary", "metadata": {}, "outputs": [], "source": [ "# Note that the structure has no hydrogens. This will interfere with quantum calcs\n", "# and possibly even with conformer generation (though it didn't this time).\n", "# Hydrogens need to be added explicitly to a structure in RDKit.\n", "m1h = Chem.AddHs(m1)\n", "m1h" ] }, { "cell_type": "code", "execution_count": null, "id": "owned-viking", "metadata": {}, "outputs": [], "source": [ "# The 2D structure doesn't make sense but that's because we need to re-run the conformer generation\n", "# to place the hydrogen atoms in the right places. In the future, you could run the \"Chem.AddHs()\"\n", "# immediately after creating the molecule from SMILES.\n", "AllChem.EmbedMolecule(m1h)\n", "m1h" ] }, { "cell_type": "code", "execution_count": null, "id": "adult-maximum", "metadata": {}, "outputs": [], "source": [ "# Ah, much better\n", "print(Chem.MolToXYZBlock(m1h))\n", "nglview.show_rdkit(m1h)" ] }, { "cell_type": "markdown", "id": "5731a70b", "metadata": {}, "source": [ "The power of cheminformatic tools becomes apparent when you want to work with big batches of molecules. For example we can create RDKit molecule objects from each molecule in the file:" ] }, { "cell_type": "code", "execution_count": null, "id": "e7896d0f", "metadata": {}, "outputs": [], "source": [ "WHO_Molecules = {}\n", "with open('WHO_EML_SMILES.csv') as csv_file:\n", " reader = csv.reader(csv_file)\n", " for line in reader:\n", " key = line[0]\n", " smi = line[1]\n", " WHO_Molecules[key] = Chem.MolFromSmiles(smi)" ] }, { "cell_type": "markdown", "id": "47bc73cb", "metadata": {}, "source": [ "Running the above command produced a few warnings for me, which is to be expected when working with large datasets, but there were no errors." ] }, { "cell_type": "markdown", "id": "bb6a50fa", "metadata": {}, "source": [ "Now that we have a list of Molecule objects, we could calculate some properties. For example, let's calculate the molecular weight of all the compounds and make a histogram. The correct API call to use can be obtained by Google searching \"rdkit molecular weight\". It gives us a function `rdkit.Chem.Descriptors.ExactMolWt` which takes a molecule object as input." ] }, { "cell_type": "code", "execution_count": null, "id": "necessary-deputy", "metadata": {}, "outputs": [], "source": [ "# For reasons unknown, we had to important ExactMolWt explicitly.\n", "from rdkit.Chem.Descriptors import ExactMolWt\n", "MolWts = []\n", "for molname, molobj in WHO_Molecules.items():\n", " MolWts.append(ExactMolWt(molobj))\n", "plt.hist(MolWts)" ] }, { "cell_type": "markdown", "id": "f38d6cd1", "metadata": {}, "source": [ "The above histogram was somewhat useful but we didn't get much info on the distribution of molecular weights other than that they are mostly under 1000. Apparently the list also contains just a few molecules with even larger weights. How about we plot a histogram with just the weights under 1000?" ] }, { "cell_type": "code", "execution_count": null, "id": "thorough-reform", "metadata": {}, "outputs": [], "source": [ "MolWts = np.array(MolWts)\n", "plt.hist(MolWts[np.where(MolWts<1000)])" ] }, { "cell_type": "markdown", "id": "ff621b91", "metadata": {}, "source": [ "Great, that's a bit more informative. We could also count the number of carbon atoms in each molecule and see if it's correlated with the molecular weight. After a bit of experimentation I figured out how to get the atomic symbols for an example molecule:" ] }, { "cell_type": "code", "execution_count": null, "id": "physical-delta", "metadata": {}, "outputs": [], "source": [ "mol_a = WHO_Molecules['abacavir']\n", "for atom in mol_a.GetAtoms():\n", " print(atom.GetSymbol())" ] }, { "cell_type": "code", "execution_count": null, "id": "81d21e78", "metadata": {}, "outputs": [], "source": [ "# In Python 3.6 onwards (we are using 3.7 I think), dictionary keys are ordered.\n", "# In past versions and possibly other Python implementations, the keys may not be ordered.\n", "# Therefore, it is a good idea to make a sorted list of keys before making element-wise\n", "# comparisons between properties.\n", "NumCarbons = []\n", "MolWts = []\n", "MolNames = sorted(list(WHO_Molecules.keys()))\n", "for molname in MolNames:\n", " molobj = WHO_Molecules[molname]\n", " MolWts.append(ExactMolWt(molobj))\n", " num_C = 0\n", " for atom in molobj.GetAtoms():\n", " if atom.GetSymbol() == 'C':\n", " num_C += 1\n", " NumCarbons.append(num_C)\n", "plt.scatter(MolWts, NumCarbons)\n", "# I set the plot limits to focus on the MW < 1000 molecules.\n", "plt.xlim(0, 1000)\n", "plt.ylim(-2, 50)" ] }, { "cell_type": "markdown", "id": "3155bb29-a054-4f23-bf79-3ec35ef24bb5", "metadata": {}, "source": [ "One of the nifty features of RDKit is that it can draw multiple molecules in a grid image. By specifying the legends keyword argument, we can put a text label beneath each molecule. The WHO essential medicine list has too many molecules to draw in a single image, so the following command draws only a few dozen." ] }, { "cell_type": "code", "execution_count": null, "id": "2eae1790-c09c-43ca-8695-20b08d83df9c", "metadata": {}, "outputs": [], "source": [ "Chem.Draw.MolsToGridImage(list(WHO_Molecules.values()), legends=list(WHO_Molecules.keys()))" ] }, { "cell_type": "markdown", "id": "06a9a9e5-879a-481e-a889-dc2021a27ba7", "metadata": {}, "source": [ "One of the most powerful tools in cheminformatics is the ability to \"search\" molecules for chemical patterns. The SMARTS language is the language for specifying a chemical pattern. The molecule represented by a SMILES string is usually, but not always, matched by the same string when used as a SMARTS. Here we will create a RDKit molecule from a SMARTS string for a benzene ring, and apply it as a pattern to a molecule." ] }, { "cell_type": "code", "execution_count": null, "id": "0830ff52-3af0-43a8-b5fe-86fffb5cfa6c", "metadata": {}, "outputs": [], "source": [ "benz = Chem.MolFromSmarts('c1ccccc1')\n", "benz" ] }, { "cell_type": "markdown", "id": "05ec99e8-c949-4d05-bc99-a111e0644bff", "metadata": {}, "source": [ "The depiction of the SMARTS benzene is slightly different from the SMILES. Now we can call the `mol.GetSubstructMatches()` method, passing the benzene SMARTS molecule as the argument, and it will return the atom indices for the matching atom groups." ] }, { "cell_type": "code", "execution_count": null, "id": "8e8db0f3-f300-4535-8159-1c69c67f9a67", "metadata": {}, "outputs": [], "source": [ "WHO_Molecules['amodiaquine'].GetSubstructMatches(benz)" ] }, { "cell_type": "markdown", "id": "d4c47627-a88c-425d-886a-2165b81220f8", "metadata": {}, "source": [ "RDKit also conveniently highlights the portions of the molecule corresponding to a positive match." ] }, { "cell_type": "code", "execution_count": null, "id": "f29bedca-74ec-4eec-aed5-d90959ca477c", "metadata": {}, "outputs": [], "source": [ "WHO_Molecules['amodiaquine']" ] }, { "cell_type": "markdown", "id": "606779ee-48c3-4430-95b1-44264da637be", "metadata": {}, "source": [ "Another great feature of RDKit is the ability to find common substructures in molecules. This can be very important for chemistry research because some compounds are \"relatives\" of one another in the sense that they have similar functions, or have the same molecular \"core\", or both. For example let's take the antiviral drugs abacavir and aciclovir. We can see that there is a fused ring structure corresponding to a purine that is shared by both. " ] }, { "cell_type": "code", "execution_count": null, "id": "6b6a3224-d51c-475f-9e2b-0ceb9e3bb265", "metadata": {}, "outputs": [], "source": [ "name1 = 'betamethasone'\n", "name2 = 'budesonide'\n", "mol1 = WHO_Molecules[name1]\n", "mol2 = WHO_Molecules[name2]\n", "Chem.Draw.MolsToGridImage([mol1,mol2], legends=['name1','name2'])" ] }, { "cell_type": "markdown", "id": "7476be49-20c4-4931-9f6b-0c6e32768356", "metadata": {}, "source": [ "Let's see if RDKit can find the common substructure for us. This requires importing a submodule of RDKit called `rdFMCS`, where the MCS stands for \"Maximum Common Substructure\"." ] }, { "cell_type": "code", "execution_count": null, "id": "cfb58733-79ec-4017-a906-94d52fb4fcbd", "metadata": {}, "outputs": [], "source": [ "from rdkit.Chem import rdFMCS\n", "# The return value of FindMCS() is a MCSResult object.\n", "# The MCSResult can be converted to a SMARTS string, which could then be converted to a molecule object.\n", "# I'm not sure why RDKit doesn't simply return the molecule object directly, but it's not so bad.\n", "mcs = rdFMCS.FindMCS([mol1,mol2])\n", "mcsMol = Chem.MolFromSmarts(mcs.smartsString)\n", "mcsMol" ] }, { "cell_type": "markdown", "id": "b96f33c0-7b90-4386-a00e-d3a0bc4c77f9", "metadata": {}, "source": [ "Now that we have the maximum common structure as a Molecule object, we can draw the two molecules together with the common structure highlighted." ] }, { "cell_type": "code", "execution_count": null, "id": "63877cba-6b0f-454d-8a6e-07634c30da2e", "metadata": {}, "outputs": [], "source": [ "mol1_matches = mol1.GetSubstructMatch(mcsMol)\n", "mol2_matches = mol2.GetSubstructMatch(mcsMol)\n", "Chem.Draw.MolsToGridImage([mol1, mol2], legends=[name1,name2], highlightAtomLists=[mol1_matches, mol2_matches])" ] }, { "cell_type": "code", "execution_count": null, "id": "f55cd6f8-71bc-4399-9ed4-e5e2124ae7d6", "metadata": {}, "outputs": [], "source": [ "Match_Molecules = {}\n", "for compound_name, mol in WHO_Molecules.items():\n", " if len(mol.GetSubstructMatch(mcsMol)) > 0:\n", " Match_Molecules[compound_name] = mol\n", "\n", "Chem.Draw.MolsToGridImage(list(Match_Molecules.values()), legends=list(Match_Molecules.keys()))" ] }, { "cell_type": "markdown", "id": "e28592e6-8c4d-44a4-9c93-f89b7fc3be96", "metadata": {}, "source": [ "With these methods in hand, you should now be able to compare pairwise similarities between long lists of molecules, or mine a large dataset of molecules for those that contain functional groups of interest. The possibilities are endless!" ] }, { "cell_type": "markdown", "id": "b08435bd", "metadata": {}, "source": [ "#### So far, we have only worked on the SMILES strings. Now let's see how we can incorporate some experimental data. \n", "\n", "There are lots of kinds of experimental data, and here we will show how to get some crystal structures into the notebook. This will require doing some work outside of the notebook as well.\n", "\n", "First, pick a SMILES string of interest. We will look at clomipramine, whose SMILES string is CN(C)CCCN1c2ccccc2CCc3ccc(Cl)cc13 . You can get the SMILES strings from the .csv file either by opening it directly as a spreadsheet in Excel or other software, or by printing the keys of the `WHO_Molecules` dictionary from above.\n", "\n", "Connect to \"CCDC Access Structures\" at the link: https://www.ccdc.cam.ac.uk/structures/ . Make sure you are on campus or using the [library VPN](https://www.library.ucdavis.edu/service/connect-from-off-campus/), otherwise you will hit a paywall at the next step. \n", "\n", "![Structure search](https://i.imgur.com/dkbIz0t.png)\n", "\n", "Click on the \"Structure Search\" tab. Click the \"Advanced ↓\" button at the lower left. Click the \"Substructure\" radio button next to \"Match condition:\", and paste the SMILES string for your molecule into the box that says SMARTS. (SMARTS is a language for searching for molecular structures, and generally a SMILES string is a valid SMARTS string). You might get an error for SMILES strings that contain \"@\" characters that specify stereochemistry. In these cases, remove the \"@\" signs from your string using a text editor, or a command line tool such as sed.\n", "\n", "![Search results](https://i.imgur.com/q1k4J2Z.png)\n", "\n", "After the search, you should see several results pop up. Crystal structues often contain more than just the molecule of interest. For example, the first entry that popped up in my search, which had the identifier BUXKIR, had the molecule of interest inside of a cyclodextrin ring. Some molecules may be similar but not identical to the molecule you're searching for. Look for a structure where the details include the name of the drug. For example, CIMPRA lists \"Chlorimipramine hydrochloride\" under Synonyms, and a quick Wikipedia search shows that Chlorimipramine and clomipramine are synonymous. That should be close enough for us. Download the CIF file to your hard drive by clicking Download -> Download Current Entry. Move the CIF file to the same folder where you have your notebook. \n", "\n", "Before you can load a structure into RDKit, another cheminformatic tool called OpenBabel is needed to convert the .cif file into a file format that RDKit can read (here we will use .mol2). In the folder where you have your notebook and the .mol2 file, run OpenBabel in your terminal (not in the notebook!) as follows:\n", "\n", "```\n", "(che155) $ obabel -icif 1125716.cif -omol2 -O 1125716.mol2\n", "1 molecule converted\n", "```\n", "\n", "Now you should be able to load the .mol2 file into your notebook. The default behavior is to remove hydrogen atoms from the structure, so we pass in the keyword argument to keep hydrogens ([see RDKit documentation](https://www.rdkit.org/docs/source/rdkit.Chem.rdmolfiles.html#rdkit.Chem.rdmolfiles.MolFromMol2File).)" ] }, { "cell_type": "code", "execution_count": null, "id": "9ebbe8f7", "metadata": {}, "outputs": [], "source": [ "mol_exp = Chem.MolFromMol2File('1125716.mol2', removeHs=False)" ] }, { "cell_type": "markdown", "id": "f2f61787", "metadata": {}, "source": [ "We can also use NGLView to show the 3D structure of the molecule here:" ] }, { "cell_type": "code", "execution_count": null, "id": "4e8c2062", "metadata": {}, "outputs": [], "source": [ "nglview.show_rdkit(mol_exp)" ] }, { "cell_type": "markdown", "id": "f920d8d5", "metadata": {}, "source": [ "An interesting question to ask is, how well does the computationally generated conformation agree with the experimental one? This touches on the cutting-edge field of crystal structure prediction, which is very important to industry as well as basic science.\n", "\n", "To look at this, we can create another RDKit molecule object directly from the SMILES string with no additional experimental data, generate a number of conformers, and compare the results to the above. The API call `AllChem.EmbedMultipleConfs` is how we ask RDKit to generate multiple conformers for a molecule ([see RDKit documentation](https://www.rdkit.org/docs/source/rdkit.Chem.rdDistGeom.html#rdkit.Chem.rdDistGeom.EmbedMultipleConfs).)" ] }, { "cell_type": "code", "execution_count": null, "id": "40e9cde4", "metadata": {}, "outputs": [], "source": [ "mol_smi = Chem.MolFromSmiles('CN(C)CCCN1c2ccccc2CCc3ccc(Cl)cc13')\n", "mol_smi_h = AllChem.AddHs(mol_smi)\n", "# The random seed argument is added so that multiple runs will give the same result.\n", "AllChem.EmbedMultipleConfs(mol_smi_h, 10, randomSeed=8282)" ] }, { "cell_type": "markdown", "id": "4ea45714", "metadata": {}, "source": [ "Using NGLView we can visualize each 3D conformation. The keyword argument for controlling which conformer to show is `conf_id` (I figured this out using `help(nglview.show_rdkit)`). I don't know whether it's possible change the conformer being viewed with a slider." ] }, { "cell_type": "code", "execution_count": null, "id": "73a92c8c", "metadata": { "scrolled": true }, "outputs": [], "source": [ "nglview.show_rdkit(mol_smi_h, conf_id=3)" ] }, { "cell_type": "markdown", "id": "ccc225e2", "metadata": {}, "source": [ "A powerful feature of RDKit is the ability to perform structure alignment. This is a calculation in which the structures of two conformations are compared by applying a translation and rotation to one confomation, called the probe, in order to minimize the distances to a reference conformation. After applying the rotation and translation, the remaining Euclidean distances between the atoms of the probe and reference structure are a measure of the \"internal\" structural differences due to differences in bond lengths, etc. By taking the RMS of all the displacements for each atom, one can calculate the \"RMSD\" which is a single number that measures the difference between two structures.\n", "\n", "In order to rotate two structures to minimize the distances between corresponding atoms, both structures need to have the same atom ordering, which is actually not true for our experimental structure and the generated one. What's more, the atoms of the probe and reference structure need to be the same, and this is not strictly true for our example either, because our experimental structure is actually a hydrochloride salt, containing an extra HCl that the generated structure doesn't have. RDKit takes care of both of these problems by finding a substructure that the probe and reference structures have in common; the substructure matching also generates a mapping of the atom orderings between one structure and the other, so we are able to get a result without worryinng." ] }, { "cell_type": "code", "execution_count": null, "id": "71e4c28e", "metadata": {}, "outputs": [], "source": [ "for i in range(mol_smi_h.GetNumConformers()):\n", " print(Chem.rdMolAlign.AlignMol(mol_smi_h, mol_exp, i))" ] }, { "cell_type": "markdown", "id": "42940a97", "metadata": {}, "source": [ "The interpretation of RMSD depends on the situation, but here we use a very rough rule of thumb that a RMSD value of 1 Angstrom or less indicates a high level of agreement for molecular structures. Above we see that conformer number 6 has the best level of agreement with experiment. Let's plot them together to look at how close the structures are.\n", "\n", "Using the `update_representation()` method, we are able to change the color of each structure being shown (called the \"component\"). The structure shown in transparent red is the experimental structure." ] }, { "cell_type": "code", "execution_count": null, "id": "1a3e684c", "metadata": {}, "outputs": [], "source": [ "view = nglview.show_rdkit(mol_smi_h, conf_id=5)\n", "view.add_component(mol_exp)\n", "view.update_representation(component=1, color='red', opacity=0.5)\n", "view" ] }, { "cell_type": "code", "execution_count": null, "id": "96223986-bb1a-48b6-87ab-374ae79b351f", "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.9.13" } }, "nbformat": 4, "nbformat_minor": 5 }