{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I have previously talked/posted about ways to construct RDKit molecules more quickly. This post revisits that topic.\n",
    "\n",
    "By default the RDKit does a lot of work when constructing a molecule. The idea here is to set things up so that we only have to do that work once for a set of molecules that we're going to work with repeatedly.\n",
    "\n",
    "There's also a diversion into thinking about what chemistry information is actually needed for things like substructure searching and tuning the molecule construction to only perceive that information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'2016.09.1.dev1'"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from rdkit import Chem\n",
    "from rdkit import RDConfig\n",
    "import os,gzip\n",
    "from rdkit import rdBase\n",
    "rdBase.rdkitVersion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "Start by reading in a set of ChEMBL molecules that we've used before and then reducing that to 50K examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "234681"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ind = [x.strip().split() for x in open('../data/chembl16_2010-2012.smi')]\n",
    "len(ind)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import random\n",
    "random.seed(0xf00d)\n",
    "random.shuffle(ind)\n",
    "ind = ind[:50000]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ms = [Chem.MolFromSmiles(x[0]) for x in ind]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "49994"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ms = [x for x in ms if x is not None]\n",
    "len(ms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's get RDKit-generated representations of the molecules:\n",
    "- A molecule's `ToBinary()` method returns a serialized form of the molecule that can be saved and then later efficiently converted back into a molecule.\n",
    "- RDKit SMILES"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pkls = [x.ToBinary() for x in ms]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "smis = [Chem.MolToSmiles(x,True) for x in ms]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Timing of standard parsing\n",
    "\n",
    "How long does it take to generate the molecules from SMILES?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 3: 15.5 s per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit [Chem.MolFromSmiles(x) for x in smis]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And from the binary format?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 3: 1.96 s per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit [Chem.Mol(x) for x in pkls]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That dramatic difference isn't really that surprising given that the format was designed to be easy to parse and that no chemistry perception needs to be done.\n",
    "\n",
    "So the binary format is really efficient, but it's not human readable and it's only useable with the RDKit. It would be cool to be able to take advantage of the portability and readability of SMILES but to have the processing not take quite so long."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an aside: another benefit of using SMILES, for at least some applications, is that they are an order of magnitude smaller than the binary format:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "binary: 22456727\n",
      "smiles: 3029371\n"
     ]
    }
   ],
   "source": [
    "print(\"binary:\",sum(len(x) for x in pkls))\n",
    "print(\"smiles:\",sum(len(x) for x in smis))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we turn off the chemistry perception when parsing from the SMILES, we can see what a difference it makes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 3: 2.47 s per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit [Chem.MolFromSmiles(x,sanitize=False) for x in smis]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Unfortunately, with no sanitization at all done, these molecules aren't that useful in many of the RDKit algorithms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Partial sanitization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the standard operations done on reading molecules from SMILES is a call to the stereochemistry assignment code, which also removes redundant or extraneous stereochemistry specifications. This can be computationally expensive and most likely not needed when reading from an RDKit-generated canonical SMILES since that already has had the extranous and redundant specifications removed.\n",
    "\n",
    "Let's see how long it takes if we skip that part (which is part of the call to `MolFromSmiles()`) and just do a normal sanitization post-SMILES parsing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def partialSanit1(smi):\n",
    "    m = Chem.MolFromSmiles(smi,sanitize=False)\n",
    "    Chem.SanitizeMol(m,sanitizeOps=Chem.SANITIZE_ALL^Chem.SANITIZE_CLEANUP^Chem.SANITIZE_CLEANUPCHIRALITY)\n",
    "    return m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 3: 9.91 s per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit [partialSanit1(x) for x in smis]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's a solid reduction over the 15.5 seconds originally required, but is still a long way off from the 2.5 seconds for the completely unsanitized version.\n",
    "\n",
    "Since the RDKit SMILES also contains information about aromaticity, we can also skip the kekulization and aromatization steps of the sanitization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def partialSanit2(smi):\n",
    "    m = Chem.MolFromSmiles(smi,sanitize=False)\n",
    "    Chem.SanitizeMol(m,sanitizeOps=Chem.SANITIZE_ALL^Chem.SANITIZE_KEKULIZE^\\\n",
    "                     Chem.SANITIZE_SETAROMATICITY^Chem.SANITIZE_CLEANUP^Chem.SANITIZE_CLEANUPCHIRALITY)\n",
    "    return m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 3: 8.04 s per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit [partialSanit2(x) for x in smis]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even better.\n",
    "\n",
    "We are still calling the ring-finding code, and sometimes we don't need information about rings (for example, all substructure queries from SMILES and many queries from SMILES), so what if we skip that?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def partialSanit3(smi):\n",
    "    m = Chem.MolFromSmiles(smi,sanitize=False)\n",
    "    Chem.SanitizeMol(m,sanitizeOps=Chem.SANITIZE_ALL^Chem.SANITIZE_KEKULIZE^\\\n",
    "                     Chem.SANITIZE_SETAROMATICITY^Chem.SANITIZE_CLEANUP^Chem.SANITIZE_CLEANUPCHIRALITY^\\\n",
    "                    Chem.SANITIZE_SYMMRINGS)\n",
    "    return m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 3: 3.87 s per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit [partialSanit3(x) for x in smis]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we're just concerned about doing substructure searches or generating RDKit fingerprints, this is still doing some extra work. Let's go to the bare minimum of sanitization: only updating the explicit and implicit valences of the atoms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def partialSanit4(smi):\n",
    "    m = Chem.MolFromSmiles(smi,sanitize=False)\n",
    "    m.UpdatePropertyCache()\n",
    "    return m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 3: 2.71 s per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit [partialSanit4(x) for x in smis]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "That's pretty fast and those molecules are actually useful for many, many calculations.\n",
    "\n",
    "We can add some ring information by calling `FastFindRings()`. This algorithm provides reliable information about whether or not atoms or bonds are in rings - and so can help with basic ring-membership queries on atoms and bonds in SMARTS or for generating Pattern fingerprints or standard Morgan fingerprints- but doesn't help with the number of smallest rings that an atom/bond is in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def partialSanit5(smi):\n",
    "    m = Chem.MolFromSmiles(smi,sanitize=False)\n",
    "    m.UpdatePropertyCache()\n",
    "    Chem.FastFindRings(m)\n",
    "    return m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 3: 3.46 s per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit [partialSanit5(x) for x in smis]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I think those last two are pretty good. By thinking about the information we need and starting from a reliable SMILES we can get molecules that are useful for many RDKit operations much more quickly than the default."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python [default]",
   "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}