RDKit_Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RDKit Tutorial, October 2020, Bender Group"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial we will cover basic functionalities of RDKit and Jupyter Notebooks.\n",
    "\n",
    "You should have followed the instructions outlined in the \"Getting Started with RDKit.pdf\" file.\n",
    "\n",
    "First, click on this \"cell\" and see that it is a \"Markdown\" cell. This means that here you can write text to explain and document your code.\n",
    "\n",
    "Now, click on the next cell and notice it says \"Code\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import packages\n",
    "# (to annotate in a \"code\" cell you must use a hash symbol)\n",
    "import pandas as pd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With Jupyter notebooks it is easy to perform exploratory analysis and test your code. For example, say you are creating a function to add 2 to an input number..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# declare new function\n",
    "def add_two(num):\n",
    "    print(num+3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# apply the function\n",
    "add_two(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hmm, looks like there's something wrong here. Instead of having to diagnose the issue, attempt to fix it and run the entire script again, we can simply modify and run just a couple of cells, separately from the rest of the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# declare new function\n",
    "def add_two(num):\n",
    "    print(num+2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# apply the function\n",
    "add_two(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That was easy!\n",
    "\n",
    "With Jupyter notebooks you can also display dataframes to ensure that everything is as it should be"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a df\n",
    "df = pd.DataFrame(np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]),\n",
    "                   columns=['a', 'b', 'c'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# display df\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# display first two rows of df\n",
    "df.head(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also display images and do a bunch of other cool stuff, but we will explore these as we go through the RDKit tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Getting Started With RDKit\n",
    "\n",
    "All material is from here https://www.rdkit.org/docs/GettingStartedInPython.html which gives a nice overview of RDKit functionalities. Thus, I give full credit to the RDKit authors and don't try to claim this tutorial as my own! (I just put it into a Jupyter notebook...)\n",
    "\n",
    "Most of the RDKit functionality is in the module rdkit.Chem. Let's import it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit import Chem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can construct molecule objects using different approaches"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = Chem.MolFromSmiles('Cc1ccccc1')\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = Chem.MolFromMolFile('input.mol') # .mol file\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stringWithMolData=open('input.mol','r').read()\n",
    "m = Chem.MolFromMolBlock(stringWithMolData)\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the import does not work, you will return None (good for error catching), e.g.:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = Chem.MolFromMolFile('invalid.mol')\n",
    "print(m is None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The error messages are often very informative:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m2 = Chem.MolFromSmiles('c1cc1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reading sets of molecules\n",
    "Reading groups of molecules is done using a Supplier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "suppl = Chem.SDMolSupplier('5ht3ligs.sdf') # sdf\n",
    "for mol in suppl:\n",
    "    if mol is None: continue # catch erroneous mols\n",
    "    print(mol.GetNumAtoms()) # print the number of atms in each mol object"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lists of molecules can also be generated:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mols = [x for x in suppl]\n",
    "len(mols)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also directly access the suppl object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "suppl[0].GetNumAtoms()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Writing molecules\n",
    "You can also convert molecule objects back to text"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = Chem.MolFromMolFile(\"input.mol\") # this provides CANONICAL smiles\n",
    "Chem.MolToSmiles(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also kekulize molecules before converting back to SMILES"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Chem.Kekulize(m)\n",
    "Chem.MolToSmiles(m,kekuleSmiles=True) # NOT canonical"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also save as MDL MOL blocks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.SetProp(\"_Name\",\"test\") # you can set a name so that it shows up in the MOL block\n",
    "print(Chem.MolToMolBlock(m))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generating coordinates\n",
    "It is possible to generate 2D and 3D coordinates too, with the 'AllChem' module"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#2D\n",
    "from rdkit.Chem import AllChem\n",
    "AllChem.Compute2DCoords(m)\n",
    "print(Chem.MolToMolBlock(m))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#3D\n",
    "m2 = Chem.AddHs(m) # add hydrogens first\n",
    "AllChem.EmbedMolecule(m2,randomSeed=0xf00d) # random seed for reproducibility\n",
    "print(Chem.MolToMolBlock(m2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove hydrogens\n",
    "m3 = Chem.RemoveHs(m2)\n",
    "print(Chem.MolToMolBlock(m3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write to file\n",
    "print(Chem.MolToMolBlock(m3),file=open('test.mol','w+'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating 3D conformations\n",
    "RDKit can also generate 3D conformations and do some energy minimisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = Chem.MolFromSmiles('C1CCC1OC') # first read in a mol object\n",
    "m2 = Chem.AddHs(m) # add Hs\n",
    "AllChem.EmbedMolecule(m2) # make 3D molecule (ETKDG method)\n",
    "res = AllChem.MMFFOptimizeMoleculeConfs(m2) #MMFF94 optimisation\n",
    "print(res) # (not_converted,energy) if not_converged=0, then the minimsisation converged"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To create multiple conformers, run ETKDG e.g. 10 times\n",
    "cids = AllChem.EmbedMultipleConfs(m2,numConfs=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now you can align the conformers and see the RMS values between them\n",
    "rmslist = []\n",
    "AllChem.AlignMolConformers(m2,RMSlist=rmslist)\n",
    "rmslist # RMS values between first conformer and all others"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can also specify RMS between two specific confomers\n",
    "rms = AllChem.GetConformerRMS(m2,1,9,prealigned=True)\n",
    "rms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Writing sets of molecules\n",
    "You can write multiple molecules, e.g. an SDF file, with the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = Chem.SDWriter('test.sdf')\n",
    "for m in mols:\n",
    "    w.write(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Drawing Molecules\n",
    "You can create images from molecules with the rdkit.Chem.Draw package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "suppl = Chem.SDMolSupplier('cdk2.sdf') # read in .sdf file\n",
    "ms = [x for x in suppl if x is not None] # make into list of mol objects\n",
    "for m in ms: tmp=AllChem.Compute2DCoords(m) # for each mol, compute the 2D coords\n",
    "from rdkit.Chem import Draw # import package\n",
    "Draw.MolToFile(ms[0],'cdk2_mol1.o.png') # draw mol file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is even possible to create an image grid out of a set of molecules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img=Draw.MolsToGridImage(ms[:8],molsPerRow=4,subImgSize=(200,200),legends=[x.GetProp(\"_Name\") for x in ms[:8]])   \n",
    "img.save('cdk2_molgrid.o.png')    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It would look even nicer if all molecules were aligned with a common core!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = Chem.MolFromSmiles('[nH]1cnc2cncnc21') # define substructure\n",
    "subms = [x for x in ms if x.HasSubstructMatch(p)] # search for mol objects which have this particular substructure\n",
    "print(len(subms)) # how many?\n",
    "AllChem.Compute2DCoords(p) # compute 2D coords of the substructure core\n",
    "for m in subms: AllChem.GenerateDepictionMatching2DStructure(m,p) # generate 2D structures aligned to this common core\n",
    "img=Draw.MolsToGridImage(subms,molsPerRow=4,subImgSize=(200,200),legends=[x.GetProp(\"_Name\") for x in subms])    \n",
    "img.save('cdk2_molgrid.aligned.o.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take a look and see the difference in the two pictures!\n",
    "Not aligned:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src = \"cdk2_molgrid.o.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aligned:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src = 'cdk2_molgrid.aligned.o.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Substructure Searching\n",
    "Substructure matching is done using query molecules from SMARTS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = Chem.MolFromSmiles('c1ccccc1O') # create mol object\n",
    "patt = Chem.MolFromSmarts('ccO') # create SMARTS query\n",
    "m.HasSubstructMatch(patt) # TRUE means it contains the substructure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can also get the atom indices for the substruct:\n",
    "m.GetSubstructMatch(patt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# display the match\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Stereochemistry is NOT used by default in substructure searches"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = Chem.MolFromSmiles('CC[C@H](F)Cl')\n",
    "print(m.HasSubstructMatch(Chem.MolFromSmiles('C[C@H](F)Cl')))\n",
    "print(m.HasSubstructMatch(Chem.MolFromSmiles('C[C@@H](F)Cl')))\n",
    "print(m.HasSubstructMatch(Chem.MolFromSmiles('CC(F)Cl')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Change this behaviour with the useChirality argument"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.HasSubstructMatch(Chem.MolFromSmiles('C[C@H](F)Cl'),useChirality=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.HasSubstructMatch(Chem.MolFromSmiles('C[C@@H](F)Cl'),useChirality=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Structure transformations\n",
    "Simple modifications of compounds can also be carried out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = Chem.MolFromSmiles('CC(=O)O')\n",
    "patt = Chem.MolFromSmarts('C(=O)[OH]')\n",
    "rm = AllChem.DeleteSubstructs(m,patt) # DELETE substructure\n",
    "Chem.MolToSmiles(rm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m # original"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "patt # sub"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rm # removed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "repl = Chem.MolFromSmiles('OC') # replacement\n",
    "patt = Chem.MolFromSmarts('[$(NC(=O))]') # pattern\n",
    "m = Chem.MolFromSmiles('CC(=O)N') # query\n",
    "rms = AllChem.ReplaceSubstructs(m,patt,repl) # replace patt in m with repl\n",
    "Chem.MolToSmiles(rms[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m # original"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rms[0] # replaced "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate Murcko Scaffolds\n",
    "Quickly get Murcko scaffolds with RDKit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit.Chem.Scaffolds import MurckoScaffold # import package\n",
    "cdk2mols = Chem.SDMolSupplier('cdk2.sdf') # load SDF file\n",
    "m1 = cdk2mols[0] # pick out first mol\n",
    "core = MurckoScaffold.GetScaffoldForMol(m1) # get Murcko Scaffold\n",
    "Chem.MolToSmiles(core) # convert to smiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "core"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate Morgan Fingerprints (ECFP4) and Calculate Similarity\n",
    "Morgan fingerprints can also be generated, useful for calculating similarity between molecules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit.Chem import AllChem # AllChem is needed\n",
    "from rdkit import DataStructs # for similarity \n",
    "m1 = Chem.MolFromSmiles('Cc1ccccc1')\n",
    "fp1 = AllChem.GetMorganFingerprint(m1,2) # radius 2\n",
    "fp1 # by default, it is a sparse integer vector object"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also generate the fingerprints as bit vectors, which is probably what you want:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fp1 = AllChem.GetMorganFingerprintAsBitVect(m1,2,nBits=1024) # 1024-length ECFP4 fingerprints\n",
    "fp1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These fingerprints can be used to compute Tanimoto similarities!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m2 = Chem.MolFromSmiles('Cc1ncccc1')\n",
    "fp2 = AllChem.GetMorganFingerprintAsBitVect(m2,2,nBits=1024)\n",
    "DataStructs.TanimotoSimilarity(fp1,fp2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the fingerprints as 0 and 1 for modelling, do the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fp1.ToBitString()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fp2.ToBitString()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also visualise the fingerprints by highlighting the atom environments which define a particular bit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit.Chem import Draw # Draw package\n",
    "mol = Chem.MolFromSmiles('c1ccccc1CC1CC1') # get mol object\n",
    "bi = {} # initialise empty dict\n",
    "fp = AllChem.GetMorganFingerprintAsBitVect(mol, radius=2, bitInfo=bi) # when generating fp, put bit info into dict\n",
    "print(bi[872]) # get a particular bit, e.g. bit number 872\n",
    "mfp2_svg = Draw.DrawMorganBit(mol, 872, bi, useSVG=True) # use this function to draw the bit\n",
    "mfp2_svg # blue = central atom, yellow = aromatic atoms, grey = aliphatic ring atoms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mol # see the full mol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Picking Diverse Molecules Using Fingerprints\n",
    "\n",
    "A common task is to pick a small subset of diverse molecules from a larger set. The RDKit provides a number of approaches for doing this in the rdkit.SimDivFilters module. The most efficient of these uses the MaxMin algorithm. \n",
    "\n",
    "\n",
    "1. Generate descriptors (e.g. Morgan fingerprints) for all the molecules, both any initial seeds plus those to pick from (the candidate pool). For large sets this will be slow, but cannot be avoided.\n",
    "2. If there are no initial seeds select an initial molecule at random from the candidate pool and it becomes the sole member of the picked set.\n",
    "3. From the molecules in the candidate pool find the one that has the maximum value for its minimum distance to molecules in the picked set (hence the MaxMin name), calculating and recording the distances as required. This molecule is the most distant one to those already picked so is transferred to the picked set.\n",
    "4. Iterate back to step 3 until your are complete (e.g. you have picked the required number of molecules).\n",
    "\n",
    "Start by reading in a set of molecules and generating Morgan fingerprints:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import\n",
    "from rdkit.Chem.rdMolDescriptors import GetMorganFingerprint\n",
    "from rdkit import DataStructs\n",
    "from rdkit.SimDivFilters.rdSimDivPickers import MaxMinPicker\n",
    "\n",
    "# load .sdf file\n",
    "ms = [x for x in Chem.SDMolSupplier('actives_5ht3.sdf')]\n",
    "\n",
    "# remove None\n",
    "while ms.count(None): ms.remove(None)\n",
    "    \n",
    "# generate fingerprints\n",
    "fps = [GetMorganFingerprint(x,2) for x in ms] # ECFP4\n",
    "nfps = len(fps)\n",
    "nfps # 180 compounds in total"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The algorithm requires a function to calculate distances between objects, we’ll do that using TanimotoSimilarity:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def distij(i,j,fps=fps):\n",
    "    return 1-DataStructs.TanimotoSimilarity(fps[i],fps[j]) # 1 - Similarity = Distance!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now create a picker and grab a set of 10 diverse molecules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "picker = MaxMinPicker() # use the MaxMin picker\n",
    "pickIndices = picker.LazyPick(distij,nfps,10,seed=23)#we are telling it to use Tc, pick 10 diverse molecules\n",
    "list(pickIndices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This gives us the Indices, but to get the molecules themselves we do this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "picks = [ms[x] for x in pickIndices]\n",
    "picks # now we have 10 mol objects of diverse molecules!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate Descriptors\n",
    "With RDKit you can also calculate a myriad of descriptors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rdkit.Chem import Descriptors\n",
    "m = Chem.MolFromSmiles('c1ccccc1C(=O)O')\n",
    "# TPSA\n",
    "Descriptors.TPSA(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#LogP\n",
    "Descriptors.MolLogP(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "List of all available descriptors can be found here: https://www.rdkit.org/docs/GettingStartedInPython.html#list-of-available-descriptors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fin."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "rdkit",
   "language": "python",
   "name": "rdkit"
  },
  "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.6.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}