Notebooks/Stiefl_RDKitPh4FullPublication.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3D pharmacophores in the rdkit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The RDKit features functionality to align molecules to pharmacophoric features in 3D. This can be used for a simple alignment as well as for a fast pharmacophore searching. Pharmcacophore features are based on so-called feature factories ( a sample file \"BaseFeatures.fdef\" can be found in the Data directory ). \n",
    "\n",
    "The method is a three step procedure:\n",
    "1. conversion of the 3D pharmacophore to a boundary matrix\n",
    "2. matching of the pharmacophore to the boundary matrix of the target molecule using triangle smoothing\n",
    "3. embedding of the molecule to the pharmacophore using distance geometry (UFF force field)\n",
    "\n",
    "Here, we provide some guidance on how to use the 3D pharmacophore tools in the RDKit based on an example. We will create a pharmacophore from an ALK5 kinase crystal structure (PDB:1PY5) and will show how to align a sample structures extracted from the crystal structure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the crystal structure, the main interaction patterns that can be identified are the hydrogen bonds between the ligand and H283 as well as the water W5 in the backpocket. In both cases the ligand interacts as an H-bond acceptor. In addition, there is a weak H-bond donor functionality interacting with D351. Apart from the directional interactions the ligand has aromatic interactions with its chinoline as well as the backpocket pyridine.\n",
    "\n",
    "![ALK5 1PY5 binding mode](images/kinaseOverview.png)\n",
    "\n",
    "Based on that information we can now start setting up a pharmacophore (ph4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "from rdkit import RDConfig, Chem, Geometry, DistanceGeometry\n",
    "from rdkit.Chem import ChemicalFeatures, rdDistGeom, Draw, rdMolTransforms\n",
    "from rdkit.Chem.Pharm3D import Pharmacophore, EmbedLib\n",
    "from rdkit.Chem.Draw import IPythonConsole, DrawingOptions\n",
    "from rdkit.Numerics import rdAlignment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting up a pharmacophore"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the RDkit a pharmacophore is defined by an object with a set of free chemical features (like acceptors, donors etc.) and a boundary matrix that describes the location of the features with respect to each other. One can define such a pharmacophore object manually."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import os.path\n",
    "fdef = os.path.join(RDConfig.RDDataDir,'BaseFeatures.fdef')\n",
    "featFactory = ChemicalFeatures.BuildFeatureFactory(fdef)\n",
    "samplePh4Feats = [ChemicalFeatures.FreeChemicalFeature('Acceptor', Geometry.Point3D(3.877, 7.014, 1.448)),\n",
    "            ChemicalFeatures.FreeChemicalFeature('Acceptor', Geometry.Point3D(7.220, 11.077, 5.625)),\n",
    "            ChemicalFeatures.FreeChemicalFeature('Acceptor', Geometry.Point3D(4.778, 8.432, 7.805))]\n",
    "samplePcophore = Pharmacophore.Pharmacophore(samplePh4Feats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.        ,  6.71795706,  6.57525467],\n",
       "       [ 6.71795706,  0.        ,  4.20853763],\n",
       "       [ 6.57525467,  4.20853763,  0.        ]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "samplePcophore._boundsMat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When applied in drug discovery projects, pharmacophore features are usually described as fuzzy locations in space. To describe this kind of fuzziness it is possible to manually change the lower and upper bounds of the respective points in the pharmacophore boundary matrix. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "samplePcophore.setLowerBound(0,1,6.0)\n",
    "samplePcophore.setUpperBound(0,1,7.0)\n",
    "samplePcophore.setLowerBound(0,2,6.0)\n",
    "samplePcophore.setUpperBound(0,2,7.0)\n",
    "samplePcophore.setLowerBound(1,2,3.5)\n",
    "samplePcophore.setUpperBound(1,2,5.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0. ,  7. ,  7. ],\n",
       "       [ 6. ,  0. ,  5. ],\n",
       "       [ 6. ,  3.5,  0. ]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "samplePcophore._boundsMat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, mostly end users will describe a three dimensional pharmacophore using a graphical interface and then write out the details to a flat file. Here, we will show an application of this using a pharmacophore (based on the above) created using MOE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "moePh43Point = \"\"\"#moe:ph4que 2014.09\n",
    "#pharmacophore 7 tag t value *\n",
    "scheme t Unified matchsize i 0 use_Hs i 1 abspos i 0 title t $ useRval i 0 comment s $\n",
    "#feature 3 expr tt color ix x r y r z r r r ebits ix gbits ix m ix\n",
    "Acc df2f2 3.877 7.014 1.448 0.3 0 400 a64cff \n",
    "Acc df2f2 7.22 11.077 5.625 0.3 0 400 a64cff \n",
    "Don f20df2 4.778 8.432 7.805 0.3 0 400 a64cff\n",
    "#endpharmacophore\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['#moe:ph4que 2014.09',\n",
       " '#pharmacophore 7 tag t value *',\n",
       " 'scheme t Unified matchsize i 0 use_Hs i 1 abspos i 0 title t $ useRval i 0 comment s $',\n",
       " '#feature 3 expr tt color ix x r y r z r r r ebits ix gbits ix m ix',\n",
       " 'Acc df2f2 3.877 7.014 1.448 0.3 0 400 a64cff ',\n",
       " 'Acc df2f2 7.22 11.077 5.625 0.3 0 400 a64cff ',\n",
       " 'Don f20df2 4.778 8.432 7.805 0.3 0 400 a64cff',\n",
       " '#endpharmacophore']"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "moePh43Point.split(\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we need to parse the MOE ph4 content and then convert it into something useful for the RDKit (a pharmacophore object)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#########################\n",
    "# parse the ph4 content #\n",
    "#########################\n",
    "def parsePh4Content(moePh4List):\n",
    "  ph4Info = {}\n",
    "  header = []\n",
    "  content = []\n",
    "  for line in moePh4List.split(\"\\n\"):\n",
    "    if line[0] == \"#\":\n",
    "      if len(header) > 0:\n",
    "        ph4Info[header[0]] = (header,content)\n",
    "        header = []\n",
    "        content = []\n",
    "      header = line.strip().split()\n",
    "    else:\n",
    "      content.extend(line.strip().split())\n",
    "  return ph4Info\n",
    "\n",
    "############################    \n",
    "# convert the feature info #\n",
    "############################\n",
    "def convertFeatureInfo(ph4InfoDict):\n",
    "  feats = []\n",
    "  radii=[]\n",
    "  # this is the conversion of the moe ph4 nomenclature to the rdkit one\n",
    "  moePh4Conv = {\"Acc\":\"Acceptor\",\"Don\":\"Donor\",\"Aro\":\"Aromatic\",\"Hyd\":\"Hydrophobe\"}  \n",
    "  header,content = ph4InfoDict['#feature']\n",
    "  noFeats = int(header[1])\n",
    "  lenSingleFeature = int((len(header)/2)-1)\n",
    "  # loop over ph4 features and populate pharmacophore object\n",
    "  for i in range(noFeats):\n",
    "    featFamily = content[(lenSingleFeature*i)]\n",
    "    (x,y,z) = (float(content[(lenSingleFeature*i)+2]),float(content[(lenSingleFeature*i)+3]),\n",
    "               float(content[(lenSingleFeature*i)+4]))\n",
    "    r = float(content[(lenSingleFeature*i)+5])\n",
    "    feats.append(ChemicalFeatures.FreeChemicalFeature(moePh4Conv[featFamily],Geometry.Point3D(x,y,z)))\n",
    "    radii.append(r)\n",
    "  return (feats,radii)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "############################\n",
    "# create the pharmacophore #\n",
    "############################\n",
    "ph4Info = parsePh4Content(moePh43Point)\n",
    "feats,radii = convertFeatureInfo(ph4Info)\n",
    "pcophore = Pharmacophore.Pharmacophore(feats)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the pharmacophore object created we need to modfy the bounds matrix to incorporate the radii of the various ph4 features. For this, we set the respective lower and upper boundaries in the bounds matrix of the pharmacophore object. For this we calculate the distance between the features and subtract (lower bounds) or add (upper bounds) the summ of the respective feature radii (see image).\n",
    "\n",
    "![MOE PH4 distance example](images/1py5Ph4DistancesOrig.png)\n",
    "\n",
    "For example, the distance between F1:Acc and F2:Acc is 6.72 Angstrom and the corresponding feature radii are 1.0 Angstrom each. So, in order to include the fuzziness of the pharmacophores the lower bounds are then set to max(featureDistance-(F1Radii+F2Radii),0). The max() is here applied in order to avoid having negative bounds values. The same is done for the upper bounds only that in this case the same of the radii is added to the distance\n",
    "\n",
    "![Fuzzy boundary matrix](images/1py5Ph4DistanceExample.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.        ,  6.71795706,  6.57525467],\n",
       "       [ 6.71795706,  0.        ,  4.20853763],\n",
       "       [ 6.57525467,  4.20853763,  0.        ]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pcophore._boundsMat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def applyRadiiToBounds(radii,pcophore):\n",
    "  for i in range(len(radii)):\n",
    "    for j in range(i+1,len(radii)):\n",
    "      sumRadii = radii[i]+radii[j]\n",
    "      pcophore.setLowerBound(i,j,max(pcophore.getLowerBound(i,j)-sumRadii,0))\n",
    "      pcophore.setUpperBound(i,j,pcophore.getUpperBound(i,j)+sumRadii)\n",
    "applyRadiiToBounds(radii,pcophore)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.        ,  7.31795706,  7.17525467],\n",
       "       [ 6.11795706,  0.        ,  4.80853763],\n",
       "       [ 5.97525467,  3.60853763,  0.        ]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pcophore._boundsMat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Checking a molecule against a pharmacophore"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The RDKit provides a simple method to check a pharmacophore against a molecule. The idea of this is to filter out quickly molecules for which an embedding is not possible. Essentially, the molecule is analyzed for its pharmacophoric features and those are checked against the features of the pharmacophore. The output is a boolean indicating whether or not the molecule can fulfill the pharmacophore purely based on its features as well as a list of all possible mappings of the pharmacophore to the molecule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ligand = Chem.MolFromSmiles(\"c1ccc(-c2n[nH]cc2-c2ccnc3ccccc23)nc1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAACWCAYAAABNcIgQAAAlP0lEQVR4nO3dfVzN5/8H8Nc50qlE\nytK6IZFQzJKwGgu5GbFls5uv2+wxY2PE3H8JM6ORe3Oz2RjfjP0Qa2yV21ES+pIpFF/dOkN3dH/O\n+/fHWUenG051zvmcOu/n43Eecp2b63WS8+66Ptfnc4mIiMAYY4wZKLHQARhjjDEhcSFkjDFm0LgQ\nMsYYM2hcCBljjBk0LoSMMcYMGhdCxhhjBo0LIWOMMYPGhZAxxphB40LIGGPMoHEhZIwxZtC4EDLG\nGDNoXAgZY4wZNC6EjDHGDBoXQsY0SC4H1qwBOnYEmjQBRCKhEzHGXoQLIWMatHEjsH07EBoKlJYC\nvMkZY/pPxPsRMqY5nTsDX38N+PsLnYQxpi4eETKmQffuAZGRgLk54OICnDghdCLG2ItwIWRMg1q1\nAoYNAx49AtatAz76SOhEjLEX4ULIWB3l5wMLFwJLlz5rGzpU8Wf5Ihkx/w9jTO/xMULGakkuB378\nEVi0CMjKAmxtgZQUwMQEkEqBiROBs2eBtm0Vi2c8PYG7d4FXXxU4OGOsWlwIGauF8HBg9mwgKQmQ\nSID584EvvlAcE6yOXA4MHAjExQGHDgGDBuk2L2PsxXjihjE13LihKGJ+fooiOHo0kJCgmBatqQgC\niqnRzz8HysqAN98EvvtOZ5EZY2riQsjYczx8CHzyCeDurlgN2r07cPIkcOAA4Oys3mv4+wO//QY0\nawZMngysXq3dzIyx2uGpUcaqUVYGbNkCLF8OPH4M2NgAwcHA2LF1XwBz+bJiRalUqhglrlvHi2kY\n0wdcCBmrpLbHAWsjJQUYPBhITlYU1V27gKZN6/+6jLG6499HGftHQkLdjgPWRvv2QFQU0KkT8Mcf\nWfjkk89RVFSkmRdnjNUJF0Jm8P7+W3EcsEePuh8HrA1HR+DPPwFv74X44YdNePPNN5GXl6f5jhhj\nauFCyAyKSKR6A4CtW4EdOxSjvnXrgEuXgP79tZvjpZeA0NBv8c477+D06dPo06cPUlNTtdspY6xa\nfIyQGRSRqOqOELm5wIoVimOBrVrpNo9MJsPUqVOxc+dO2NnZ4cSJE+jWrZtuQzBm4LgQMoMiEgEt\nWypWa/r4ABs2AA4OwmYiIixbtgzLli2DlZUVjh07Bi8vL2FDMWZAeGqUGRQiIDsbuHkTaNcOmDBB\n6ESASCTC0qVLsWLFCmRnZ2PYsGGIjo4WOhZjBoNHhMxg5eUB9vaKi2fri3379iE4OBgnT55EK13P\n0zJmoLgQMoOUk6M4Qf7PPxUXyNYnxcXFkEgkABSjxYpatWqFhw8fChGLsUbLSOgAjOlSeV2xsAC8\nvYHdu4XNU53yIggojh+Wi46ORnh4uBCRGGvUeETIDFZJCWBsLHQK9fn5+WHnzp2wtbUVOgpjjQov\nlmEGRy4HXF0VJ85XlJ0NhIUB168Lk+t5zp8/DysrKy6CjGkBF0JmcMRixSkUSUmK4lfuzz+Bt98G\n9uwRKlnNgoKCEBgYKHQMxholLoTMIPXsqTiV4vLlZ22enoo/L10SJlNNzp49i7KyMri7uwsdhbFG\niQshM0jVFb2XX1acTnH5MiCT6S6L7AWdBQUFYebMmboJw5gB4kLIDFJNoz9PT+DJE8W0qS48fvwY\njo6OWLVqlUp7amoqkpOTAQCnTp3C22+/rZtAjBkgLoTMILm4AC1aVF8IAd1Nj27evBnp6emQy+Uq\n7V9++SVcXFxw8uRJ3QRhzIDx6RPMYPXvD5w+DaSnA3Z2iraICMXGudOmAZs2abf/wsJCtGvXDgUF\nBbh//z4sLS0BAJmZmXBycoKNjQ2Sk5NhZMSn+zKmTTwiZAarfPQXF/esrUcPQr9+15GT85PW+9+7\ndy+kUikCAgKURRAAtmzZguLiYkyfPp2LIGM6wCNCZrAOHcrHhg338dZbDzBr1gBle8eOHZGWlobc\n3FwYa+mMe7lcDldXV9y5cwe3b9+Gk5MTAODJkydwdHSETCZDamoqmjdvrpX+GWPP8IiQGSx394c4\ne7YrIiK+UWn39PREUVEREhIStNZ3eHg4kpKS4O/vryyCALB79248fvwYH330ERdBxnSECyEzWE5O\nTrC2tkZcxblRAD179gQAXNLiipm1a9cCAObOnatsk8lkWLt2LZo2bconzzOmQ1wImUHz8PDAw4cP\ncffuXWWb5z8HDysXSE25ePEizpw5A29vb2VfABAWFoa7d+/i3XffhYPQuwUzZkC4EDKDVl6IKo7+\nPDw8YGRkpLUR4fr16wEAs2fPVmkvHyVWbmeMaRcXQmbQyqdBK47+zMzM0LlzZyQkJODp06ca7S8l\nJQUHDx5Ex44d8dZbbynbo6OjceHCBfj4+MDDw0OjfTLGno/XZjOD1qtXLwBVjweuWrUKzZs31/iq\n0c2bN0Mmk2HmzJkQi5/9HhoSEgKAR4OMCYFPn2AGz8HBAXl5ecjJyVEpTpqWk5ODtm3bQiKR4H//\n+x/MzMwAAMnJyejUqRNcXFxw48aNKrvSM8a0i0eEzOB5eHjg2LFjuHv3Ljp06KC1fpo1a4atW7fi\n6dOnyiIIABs2bFCOErkIMqZ7PCJkBi8lJQUtW7aElZWVzvt+9OgR2rZtixYtWuDevXuQSCQ6z8CY\noePFMszgtW/fXlkEly5dqtNR2Y4dO1BQUIApU6ZwEWRMIDwiZOwfcXFxGDlyJDIzM6Gr/xZubm5I\nSUnBvXv3YGNjo5M+GWOqeETIGICioiJMmjQJ+/bt02m/0dHROHLkCBdBxgTEI0LGAMycOROOjo4I\nDAyESCQCESEsLAy5ublo06YNHBwc0KZNG5iYmAgdlTGmYVwIGQMgFourTIf6+Pjg9OnTKm0mJiaw\ns7ODra0t7Ozs0L59e7Rv3175d2dnZ1hYWDy3r19//RVz5szB3bt34eTkhG+++QZ+fn6afks6U/GY\nasXvYeVjra1atcLDhw+10pdcLkdISAi2b9+OlJQUyOVynU1vs4aPCyFjlZSPCH/77TckJibi/v37\nSE1NRVpaGlJTU5GVlVXjh6xYLMbLL7+Mtm3bwsHBAQ4ODnB0dETXrl3h6+sLALC2tsaePXswYMAA\nREZGIiAgAFKpVJdvUSvKv2/ViY6ORnh4OFasWKGVvtavX48tW7YgNDQUPXr00Or5oKzx4fMIGavB\nsGHDMGzYsGrvy87ORkpKCjIyMpCZmVnl67i4OMTExCgfP2rUKGUhtLOzA/BsdGNvb6/ldyK8r776\nCjt37sTevXuRnp5eq+fa2Nhg4sSJz33Mtm3bEBwcrLxkHmO1wYWQGaTU1FTY2tpWuwO8OpMklpaW\n8PDwqPG6oIWFhbh//z7S0tKQlpamshhmx44dGDJkCHJzc2FhYYHff/+97m+kATh//jysrKxga2uL\nb7/9FhcuXKjV83v06PHCQnjv3j1ERkZi3LhxsLOzw8aNGzF06NB6pGaGhKdGmcFJT0+Ht7c3evTo\ngQMHDlRbDLXJxcUFGzZsUE6Nzpo1C0lJSTrNoA01TY36+vrim2++gbu7Oy5duoTc3Nxqn5+dnV1t\nu6WlpXI0XVNf9vb22LFjB3x9fREZGYnJkyfXeuTJDBePCJlBkUql8PHxwf/+9z+MHj1a50UQgLIQ\nlE+N1lQYGoOzZ8+irKwM7u7uAKCy/6ImlY/+yr+nfIyQ1QoxZiDy8/PJ09OTAFBAQADJ5XJBchw5\ncoRcXFxIIpGQi4sLHT16VJAcmgKgyq2cj48PHT58WOt9PXjwgN58801q1qwZdenShSIiIjTWJ2v8\neGqUGYTi4mIMHz4cUVFRGD58OA4fPoymTZsKHQsAcPXqVRgbG8PNzU3oKHVSfqyTsYaK5w9YoyeX\nyxEQEICoqCi8/vrrOHjwoN4UwcuXL8PT0xOfffaZ0FHqzN/fHwMHDkROTo5gGXJycjB58mTMnDlT\nsAys4eJCyBq9zz//HKGhoejevTt+/fVXmJqaCh1JycPDA/3798eZM2fw888/Cx2n1kJDQ3Hq1CkQ\nEVq2bClYDjMzM5w+fRqbN2/GtWvXBMvBGiaeGmWN2jfffIO5c+eiXbt2OH/+vPIcPn0SHx+Pnj17\nws7ODomJiSp7Feqzp0+folOnTpBKpYiPj4erq6ugeQ4dOoR33nkH/fv3x8mTJwXNwhoWHhEyjTp4\n8CDc3NxgZmYGT09PnD17VrAsP/zwA+bNm4fWrVsjIiJCL4sgALz66quYOHEiUlNTsW7dOqHjqG3V\nqlVIT0/HlClTVIqgTCYTJE/5RQtOnTqFo0ePCpKBNVCCLtVhjc7o0aPp+vXrVFhYSDt37iQ7OztB\ncoSHh1PTpk3J3NycYmNjBclQGw8ePKAWLVqQubk5paenCx3nhW7dukUSiYSsra3p8ePHyvbS0lLy\n8PCgOXPmaHVVbklJCa1evZomTZqk0h4fH09NmjShDh06UFFRkdb6Z40LF0KmNXfu3KEOHTrovN8/\n//yTTE1NSSKRUFRUlM77r6uVK1cqT+3Qd/7+/gSAtm3bptK+Zs0aAkAfffSRVvuXy+XUq1cvAkAn\nTpxQue+jjz4iALR27VqtZmCNBxdCphWZmZnk4eFBYWFhRER08+ZN+vbbbyk+Pp7Kysq01u9ff/1F\nrVq1IrFYTKGhoVrrRxsKCwupXbt2JBaL6dKlS0LHqdGJEycIAHl6epJMJlO2Z2RkUPPmzally5b0\n4MEDreeIjo4mkUhErq6uVFpaqmwvH123bNmS/v77b63nYA0fF0KmcRcvXqR27drR3r17lW2rV69W\nngBtZGRErq6uNHnyZNq9ezclJCRorO/AwEACQKtXr9bYa+rSzz//TADIy8tLsBP+n6ekpIS6dOlC\nIpGIzp8/r3JfQEAAAaCQkBCd5fnggw8IAG3dulWlvXx0/dlnn+ksC2u4uBAyjdqxYwfZ2NjQH3/8\nodJ+7do1WrlyJb311lv08ssvV7k6SIcOHWjMmDG0adMmiouLU/kNvzZkMplGr2QihL59+xIA+uWX\nX4SOUsWGDRsIAI0ZM0alPSYmhsRiMXXt2rXO/3Z1kZqaSmZmZmRlZUWPHj1SthcXF1OHDh2oSZMm\ndO3aNZ3lYQ0TF0KmUZULHADKz8+v8rjMzEw6evQoBQUFkZ+fH1lZWak8p7pRY+UR0v79+6lTp05k\nampKPXv2pNOnT+vqbWrVlStXSCwWk5OTExUWFgodR0kqlZKlpSWZm5tTWlqasl0mk5GHhwcBoMjI\nSJ3n+ve//00AKDAwUKW9fHQ9cOBAnWdiDQsXQlYvRUVFlJqaWu/XKSkpoYsXL9LGjRtpzJgx5Ozs\nXKWgtm3blj744APKysoiIqJx48bRvXv3qKCggPbt20e2trb1zqEvJkyYQABo1apVQkdRmjx5MgGg\nlStXqrT/8MMPBIDefvttQXLl5+eTnZ0dNW3alBITE1Xu69evHwGg8PBwQbKxhoELIauzsrIyGj16\nNNnZ2VFSUpLGXz8/P5/OnTtH69evp9GjR5ONjQ2ZmppSSUmJyuOKiorol19+oVdeeUXjGYSSnp5O\n5ubm1Lx5c8rIyBA6DsXGxpJYLCYXFxcqLi5Wtufk5JCNjQ2ZmJhQcnKyYPl27dpFAMjPz0+lvXx0\n3blz5yo/N4yV40LI6kQul9PYsWMJAPXt25cKCgpe+Jx58+bRqFGjKDg4mM6cOUNPnjypVZ8ymYzu\n3bun0lY+WrS0tKTo6OhavZ6+W758OQGgjz/+WNAccrmcvLy8CECVnTLmzJlDAGjRokUCpVOQyWTK\nnUWOHz+ucl/56HrDhg0CpWP6jgshq5PyD8Du3btTTk6OWs957bXXqhwHfPXVV2nKlCn0ww8/0F9/\n/VWnlZLlU6POzs61fq4+e/LkCdnb21OfPn10ugClsqtXr5KRkRH5+vqqtCcmJpKxsTE5OjrS06dP\nBUr3zPnz50kkElGXLl1URn/3798nMzMzsra2puzsbOECMr3FhZDVWnBwMAGgdu3a1foqKOnp6cpF\nMr6+vmRqaqpSHM3MzMjb25s+//xz2r1793On22bMmEFSqZQKCgooNDSU2rRpU9+3pndu376tcq6e\nUG7evEm3b99WaRs5ciQBoN27dwuUqqp3332XANDmzZtV2tevX0+7du3Si+8l0z9cCFmtfPfddyQS\nicjGxqbKB2Nd5Ofn06lTp2jlypU0YsQIat26dZVFMh07dqT79+9Xee73339PdnZ2ZG5uTl5eXnTh\nwoV659FnERER5OrqSubm5jR8+HB6+PChRl+/4ve8omPHjlHnzp1JIpFQ586d6dixY0SkGIFNnTpV\nr853vHXrFrVv3155IQfG1MG7TzC1HTp0CO+99x5MTU1x8uRJeHp6aqWfzMxMxMXF4fLlyzh//jyu\nX7+OjIwMiMU1XyOeiFBQUIBmzZppJZM+aNOmDbZt2wZfX19ERkYiLCwMO3bs0Hg/IpEIFT8WrK2t\nsWfPHgwYMACRkZEICAiAVCrVeL+aIpPJ0KRJE5U2kUik/Lrie6vYXvk+ZjiMhA7AGoazZ89i7Nix\nMDIyQlhYmNaKIADY2tpixIgRGDFiBACgtLT0uUXwypUrGDFiBN5//32EhIRoLZfQynd1KP/wPn78\nOLKzs5/7HHNz83pvQly+a0d5v/b29vV6PW2rXASBZwWucuGreB8zYMINRllDceXKFWrRooXeXr/z\n6dOnZGZmRg4ODno1Tadpv/76Kzk7O5OZmRnNmjWLTExMyNjYuNqLGKh7mzx5cpV+Kn8sxMTEkIWF\nBQEgCwsLiomJ0dVb1rjK7w0AtWzZkqysrGjUqFEaOSeWNTw8ImTPdevWLQwdOhR5eXnYsmULPvjg\nA6EjVWFmZoZBgwYhLCwMcXFxWh2tCmn48OEYPnw4AODIkSM4fvw4HB0dUVZWVuNznjx5gtLS0hrv\nd3R0fGG/48aNQ2hoqHJqdPz48UhKSqr9G9BD9M9oUCqVYvXq1ZgwYQKioqIETsV0jQshq1FWVhbe\nfPNNSKVSLFy4EJ9++qnQkWrk7++PsLAwHD58uNEWQgCQy+W4evUq5s6di+nTp2P69Ola7zM3NxfA\ns2nF8r83Jq1bt0ZQUJDeT/syLRF6SMr0l4+PDwGgcePG6f2U46NHj8jIyIhcXFyEjqI1AEgsFlOH\nDh1ozZo1Gj8VANVMnRIRHTlyhFxcXEgikZCLi0uVk+obkpo+8rKzs2nBggXUt29fHSdi+oBXjbIa\nnTt3Dlu3bsWePXvqveBCFwYOHIiTJ0/ixo0bcHV1FToO0yM1LZIpb7ewsIC3tzc2b94MJycnXcdj\nAqt5KR4zKCKRSHkr17dvX4SGhjaIIggopkcB4PDhwwIn0ZwVK1YgISFB6Bg4ffo07t69K3SMOsvI\nyIBcLgcpzp1WHhss/zonJwfh4eFcBA0UF0IGACofDg3VyJEjIRKJEBYWJnQUjYiKisLixYsxYcIE\nQf9tIiIi0L9/f3zxxReCZaiP0tJSDBw4ED4+PigqKhI6DtNDXAhZo9G2bVuMGXMemZkXkJYmdJr6\nKS0txbRp0yASibBp06Zqp/Z0xdfXF97e3jh06BD++OMPwXLU1bp163Dz5k20a9cOJiYmQsdheogL\nIWtU3NxeQ1qaEf7v/4ROUj+bNm1CYmIi/vWvf8HLy0vQLCKRCGvWrIFIJEJgYOBzT9fQNxkZGVix\nYgUsLCwQHBwsdBymp7gQskblnXcUfzbkw4SZmZlYunQpzM3NsXr1amU7EWHLli3Iy8vTeaY+ffrg\nww8/xF9//YXvvvtO5/3X1cKFC5Gfn4+FCxfCxsZG6DhMXwmzWJXpq8bwI+HmRtSkCZFUKnSSupk0\naVK1O8H/9NNPBIDGjh0rSK7U1NQGtZ1RdHQ0iUQi6tSpk8pmwoxVxiNCBgAqK0Yrrx5taPz9AZkM\nOHZM6CS1Fxsbix9//BEuLi6YPXu2sj03NxezZ8+GiYkJli9fLkg2BwcHzJ49G3///Te++uorQTKo\nSy6XY9q0aSAibNq0CcbGxkJHYnqMCyEDAJVl5dTAV5D+cxZFg5seJSLMmDEDcrkca9asUfnw/vLL\nLyGVSjFnzhxBl/jPnz8fbdq0wcaNG3H79m3BcrzI7t27cfnyZfj5+WHQoEFCx2H6TsDRKNND169f\np5kzZ9LFixdV2vfv30+TJ0/W+yvMlHNyIpJIiPLyhE6ivj179hAAGjZsmEr7zZs3qWnTpuTo6EgF\nBQUCpXvmxx9/JADk7+8vdJRq5eTkUOvWrUkikdCdO3eEjsMaAC6ETMXhw4cJAE2bNk2l3d/fnwBQ\nbGysQMlqJzCQCCA6cEDoJOrJzc0lGxsbkkgkVTY89vX1JQB6s/OHTCYjT09PAkARERFCx6niiy++\nIAA0f/58oaOwBoILIVNRUFBAzZo1I3t7e5XRX/loZcGCBQKmU19cHNGXXxKlpAidRD3z5s0jADR3\n7lyV9rCwMAJAPj4+AiWr3oULF0gkElH37t2prKxM6DhKiYmJZGxsTPb29pSfny90HNZAcCFkVYwa\nNYoAqOw7l52dTcbGxtS5c2cBkzVOSUlJZGxsTHZ2dpRXYS63sLCQ2rdvT0ZGRpSQkCBgwuq99957\nBIB27twpdBSlQYMGEQD66aefhI7CGhBeLMOqqO6anS1btsQbb7yBxMRE3Lx5U6hoahGJnt0qiowE\n3NyA5s0BPz/g0SNh8lU2e/ZslJSU4Ouvv0bz5s2V7evWrUNKSgo+/vhjuLm5CZiweqtXr4aJiQkW\nLVqkF1szhYeHIyIiAv369cOYMWOEjsMaEN59glWRk5MDGxsbODo64tatW8r2b7/9Fp9++im++uor\nLFy4UMCE6hGJgIo/3W3aANu2Ab6+iqIYFgbs2CFcPgD47bffMHz4cPTt2xdnzpxRnraSnp6Ozp07\nw8TEBElJSbCyshI2aA0WLlyIr7/+GvPnz8fXX3+t1nPCwsKQmJio1mMnTJiAl19++YWPKy4uhpub\nG+7du4dLly7B3d1drddnDACvGmXVGzx4MAGg69evK9uysrJILBZTz549BUymvso/3ba2RL/+SlRc\nrPjTwUGYXBXduHGDhg4dStHR0SrtY8aMIQC0efNmgZKpJycnh6ytrWnEiBFqrygun1JV5xYXF6fW\na54+fZqMjY0pICCgPm+HGSgeEbJqbdu2DVOnTsXy5cuxePFiZbuXlxdiYmJw7949tG3bVsCEL1Z5\nRBgeDsycCWRkAFOmAFu3AoWFgsWrUWFhIfr164eCggLEx8fr/TZY6enptdrZPSYmBqmpqWo91tfX\nF5aWlmo9Njk5Gc2bN0fr1q3VzsIYAB4RsuqVj/7c3d1V2oODgwkAbdiwQaBk6nveT/fhw0Rduugs\nSpVd38tlZWXRoEGDqHnz5jRo0CDKysoiIsUpChkZGboLWA81vbf9+/dTp06dyNTUlHr27EmnT5/W\nWl/lgoKCGsVlAplu8WIZVi0bGxv06dMHV69eRUpKirJdXza/jYoCxowBiotr9zy5HLh8GZg7F5g6\nVTvZqkM1XK1n/vz56NatG1JTU9G1a1csWLAAACAWi2Fra6u7gPVQ03sLDw/H77//jkePHiEwMBAf\nfvih1voCgLi4OOwQ+qAva5CMhA7A9Je/vz8uXLiAI0eOYNasWQAAZ2dndOvWDefOncPff/8Na2tr\nnee6cgUYNUoxrTl1KvD666r3V1wtWv41keJrsRhwclI877PPdJe5Jn/88QdiY2NhYWGB2bNno3fv\n3kJH0pg9e/YAUCxkkUgksLa2RnZ2dq1fx8TEBKamps99TFFRESZNmoR9+/ZhwIABdcrLDJiQw1Gm\n31JSUggAvf766yrtS5YsIQD0/fff6zzTX38RtWpFJBYT7d+vel9hIVFams4j1Url/3LGxsZUUlJC\nREQlJSVkbGwsRCyNqO7jBP9MY1paWtLJkyfVXiRT8VbdRRwq9zVjxgwKCQmpMQdjz8OLZdhzde/e\nHQkJCUhLS1NO1cXHx8Pd3R1+fn44psMtHu7fB7y9gbQ0xUKXilObMhnw3ntAbCxw8iTQsaPOYtWK\nSCRSmdqzt7dHbGws7O3tkZ6ejt69eyMtLU3AhHVX+b2VKywsxOHDh7FkyZI6XTD8ww8/xKRJk57b\nl1gsrtI3f7QxdfHUKHsuf39/XLt2DUePHsUnn3wCQFEc33//ffj4+OgsR15eHubN+wZS6WLMmmWs\nUgTlcmDCBODQIaBfP8DBQWex6m3QoEEICQnBkiVLEBIS0qh2Spg5cyYWLVoEc3NziMVilJSUICIi\nQit9yeVy5dc1FWTGaiTcYJQ1BPHx8QSAhgwZIliGwsJCeuONNwgABQYGU+XT1aZNU6wQdXcnys0V\nJuOLoJopPyKizMxMGjhwIJmbm9PAgQMpMzNT4KS1V9N7+/7778nOzo7Mzc3Jy8uLLly4oLW+Kj+G\nsdrgqVH2Qs7Ozrh//z4ePHig9jldmlJcXAw/Pz9ERkbi3Xffxf79+9GkSRPl/cHBwLx5gIsLcO4c\n0BBPIZPJZCrviTGmW3z6BHuhESNGoHv37njw4IFO+5XL5QgICEBkZCQGDx6Mffv2qRSM/fuLsWAB\n8NJLisulNcQieOLECXTq1AkXL14UOkqdEBE+++wzxMfHCx2FsTrjESF7obKyMhgZ6f5wckJCAnr1\n6gVnZ2ecOXNGZTR66NAhfPLJFHTseATbt3uhWzedx9OIU6dOYcCAAXB3d0dcXBzE4ob1u+nu3bsx\nceJEDB06FMePHxc6DmN10rD+1zFBlBfBBw8eYPDgwWjRogUGDx6s9RFi165dERUVhRMnTqgUwZiY\nGIwfPx55eblYuvRJgy2CANC/f3+8/fbbuHr1Kvbu3St0nFrJycnB3LlzIZFIsHnzZqHjMFZnXAiZ\n2mq6Coo2vfbaa7Czs1P+/erVqxgyZAiKi4tx8OBBDB48WOsZtG3t2rWQSCRYuHAhnj59KnQctS1b\ntgxSqRQzZsxAhw4dhI7DWN0JulSHNSh2dnaU9s8Z62lpaWRvb6/T/m/dukWtW7cmkUhEu3fv1mnf\n2jZ79mwCQEFBQUJHUcu1a9fIyMioymbCjDVEfIyQqU0ikeDJkydo2rQpSktLYW5ujtjYWHTu3BkS\niUSrfWdlZeH1119HcnJyg9kPsTby8vLg4uKCvLw83Lx5E46OjkJHei5fX19ERUXhxx9/xIQJE4SO\nw1i98NQoU9tLL70EqVQKAJBKpXjppZfQo0cPNGvWDG5ubhg/fjw2bNiAy5cvq5zgrAlisRgtW7bE\nlClTGl0RBIAWLVogKCgIhYWFWLRokdBxnuvIkSOIiopCnz59MH78eKHjMFZvPCJkaps4cSJatWqF\nJUuWYPny5UhPT4elpSViYmKQkJCAsrIy5WOtra3Ru3dv9O7dG3369EGvXr3QokWLevX/9OlTmJiY\nNNpz7mQyGdzd3ZGQkIBz587B29tb6EhVlJQAfn7/RWbmVGzfvgZeXl5CR2Ks/oSdmWUNyfOuglJa\nWkoJCQm0fft2GjduHLm6upJIJFK5AoitrS35+fnRqlWr6Ny5c1RYWKh8/oEDB8jV1VW5d92ZM2eI\nqOb9+hqrqKgoAkB9+vRRe8d3XQoOVlzFZ9w4oZMwpjk8ImRak5qaipiYGMTExODixYu4cuUKCits\nCW9lZYWsrCw0bdoU7733HpYsWQJnZ2fs3bsXQUFBSE9PR0BAAKysrLBkyRIsW7YMOTk52LVrl4Dv\nSvv8/PwQHh6O//znPxrZw09T7t8HunQBjI2BpKSGeQEDxqrDhZDpTGlpKf773//i4sWLiImJQXFx\nMQ4cOFDlccnJyRgyZAju3LnTqHZnUFdSUhK6desGGxsbJCUlwczMTOhIABQbIf/nP8C6dcDMmUKn\nYUxzuBAyvZKVlQU/Pz8sWbIEI0eOrHalanFtt6VvgGbMmIGNGzdi+fLlWLx4sdBxcPYs4OMDdO2q\n2BhZgAsNMaY1vGqU6Y3Y2Fi89tprCAwMxMiRIwFUXalqbW0tZESdWbp0KVq1aoXg4GBkZGQImkUm\nU4wAiYD167kIssaHCyHTCzt37sTIkSOxY8cOjBkzRtlevl9fbm5uo9uv73ksLS2xYMECPHnyBMuX\nL1f7eTY2NhCJRCq3Pn0yIBJB5TZnjvpZdu0Crl4F3n0XGDCgDm+GMT3Hv9sxvTB58mQAULlkWn5+\nPlatWoWxY8fCwcEBvXv3bnDX46yP6dOnIyMjAzNmzFD7Oa+88gqys7NV2hwdZahwZguA2m1e3Lcv\n8M47wOrV6j+HsYaEjxEyxhgzaDw1ypieqji9qU57/fp6dqvo4EHAzQ0wMwM8PRWLZhhrbLgQMqan\niAjVTdjU1F6/vhS3yg4eBH7+GXj8GPjkE0CPTmtkTGN4apQxPScSiaotfBXbb9y4gWXLllX7fCen\nEOTkVD0oaG8PLFlS+TWrL4gAkJwMDBkC3LlTu/yM6TteLMNYIyCVSnHw4MFq7+vZcyvi4qq2d+tW\ntRDWJCsLeP99ICSkHiEZ01M8ImRMz6kzIiwsLERmZma1z3/yxBalpaZV2k1NAVfXyq9ZdUQYG6so\ngitWKK4uw1hjw4WQMT2nTiHUXF+qhXDnTmDxYuCnnwADOYWTGSAuhIzpqepWhRJRje3166tqG1H1\n7fn5gLl5vbpjTK9wIWSMMWbQ+PQJxhhjBo0LIWOMMYPGhZAxxphB40LIGGPMoHEhZIwxZtC4EDLG\nGDNoXAgZY4wZNC6EjDHGDBoXQsYYYwaNCyFjjDGDxoWQMcaYQeNCyBhjzKBxIWSMMWbQ/h+P4LXj\nHrl5mgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<rdkit.Chem.rdchem.Mol at 0x7f8ed5f00c20>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DrawingOptions.bondLineWidth=1.8\n",
    "DrawingOptions.atomLabelFontSize=14\n",
    "DrawingOptions.includeAtomNumbers=True\n",
    "ligand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "canMatch,allMatches = EmbedLib.MatchPharmacophoreToMol(ligand,featFactory,pcophore) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "canMatch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 Acceptor SingleAtomAcceptor (5,)\n",
      "0 Acceptor SingleAtomAcceptor (12,)\n",
      "0 Acceptor SingleAtomAcceptor (19,)\n",
      "1 Acceptor SingleAtomAcceptor (5,)\n",
      "1 Acceptor SingleAtomAcceptor (12,)\n",
      "1 Acceptor SingleAtomAcceptor (19,)\n",
      "2 Donor SingleAtomDonor (5,)\n",
      "2 Donor SingleAtomDonor (6,)\n"
     ]
    }
   ],
   "source": [
    "for (i,match) in enumerate(allMatches):\n",
    "  for f in match:\n",
    "    print(\"%d %s %s %s\"%(i, f.GetFamily(), f.GetType(), f.GetAtomIds()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matching a molecule to a pharmacophore"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the RDKit pharmacophore module is based on distance geometry, a molecule can be matched to a pharmacophore without actually aligning it to the pharmacophore. For this, the molecules' (smoothed) bounds matrix is updated with the distances of the pharmacophore bounds at the respective atom matches. Then, the algorithm tries to smooth the modified bounds matrix. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "boundsMat = rdDistGeom.GetMoleculeBoundsMatrix(ligand)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if canMatch:\n",
    "  failed,boundsMatMatched,matched,matchDetails = EmbedLib.MatchPharmacophore(allMatches,boundsMat,\n",
    "                                                                             pcophore,useDownsampling=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apart from an indicator if the molecule could be matched to the pharmacophore (failed), a lot of useful information is returned. The molecules' smoothed bounds matrix (boundsMatMatched) as well as information of what atoms were matched to the pharmacophore (matched)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "failed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(21, 21)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "boundsMatMatched.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Acceptor 12\n",
      "Acceptor 19\n",
      "Donor 5\n"
     ]
    }
   ],
   "source": [
    "for match in matched:\n",
    "  print(\"%s %d\"%(match.GetFamily(),match.GetAtomIds()[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computation times of bounds smoothing depend heavily on the size of the bounds matrix. To speed up things, the RDKit provides an argument to the MatchPharmacophore method (useDownsampling) which essentially reduces the bounds matrix to the relevant elements and then does the bounds smoothing on the reduced matrix. The matchDetails return value is a tuple that provides information on the size of the input bounds matrix as well as the output bounds matrix. If downsampling is not enabled both values are identical if not they differ where the second element is the reduced bounds matrix size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(21, 21)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "matchDetails"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "failed,boundsMatMatched,matched,matchDetails = EmbedLib.MatchPharmacophore(allMatches,boundsMat,\n",
    "                                                                           pcophore,useDownsampling=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "failed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(21, 13)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "matchDetails"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Embedding a molecule onto a pharmacophore ( NOT the actual alignment)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we know that the molecule should actually embed onto the pharmacophore and how (using the list of matches returned from MatchPharmacophore), embedding the molecule (with the ph4 constraints) is actually quite straightforward."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "atomMatch = [list(x.GetAtomIds()) for x in matched]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[12], [19], [5]]"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "atomMatch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To avoid really bad conformations one needs to add hydrogends at this point (it helps a bit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ligH = Chem.AddHs(ligand)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "bm,embeddings,numFail = EmbedLib.EmbedPharmacophore(ligH,atomMatch,pcophore,count=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "numFail"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(embeddings)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case none of the embeddings failed. Note - at this stage the molecule is just embedded onto the pharmacophore and not aligned to it. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aligning an embedding to a pharmacophore"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to do the the alignment code provided by the RDKit can be used. First we start by generating the coordinates of the matching pharmacophoric points - since this can be multiple points per match we cannot just use the atomIds and do the alignment based on that. For this, we create the geometric mean of the multiple atoms coordinates and use that to align to the pharmacophore features. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.3324263571321069,\n",
       " 0.33586698144527105,\n",
       " 0.32520986638642313,\n",
       " 0.32671130271598514,\n",
       " 0.3421387037292547,\n",
       " 0.3224447247298343,\n",
       " 0.3425788770183189,\n",
       " 0.3339505837944188,\n",
       " 0.31849689244806356,\n",
       " 0.3298604798363485]"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def GetTransformMatrix(alignRef,confEmbed,atomMatch):\n",
    "  alignProbe = []\n",
    "  for matchIds in atomMatch:\n",
    "    dummyPoint = Geometry.Point3D(0.0,0.0,0.0)\n",
    "    for id in matchIds:\n",
    "      dummyPoint += confEmbed.GetAtomPosition(id)\n",
    "    dummyPoint /= len(matchIds)\n",
    "    alignProbe.append(dummyPoint)\n",
    "  return (rdAlignment.GetAlignmentTransform(alignRef,alignProbe))\n",
    "\n",
    "def TransformEmbeddings(pcophore,embeddings,atomMatch):\n",
    "  alignRef = [f.GetPos() for f in pcophore.getFeatures()]\n",
    "  SSDs = []\n",
    "  for embedding in embeddings:\n",
    "    conf = embedding.GetConformer()\n",
    "    SSD,transformMatrix = GetTransformMatrix(alignRef,conf,atomMatch)\n",
    "    rdMolTransforms.TransformConformer(conf,transformMatrix)\n",
    "    SSDs.append(SSD)\n",
    "  return(SSDs)\n",
    "\n",
    "TransformEmbeddings(pcophore,embeddings,atomMatch)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All alignments seem to work pretty well (low squared deviation values). This is also true when looking at how the pharmacophoric features of the lgand fit the pharmacophoric query.\n",
    "\n",
    "![Aligned Embeddings](images/alignEmbed3Point.png)\n",
    "\n",
    "Unfortunately the query wasn't setup very well with respect to \"fixing\" the molecules. I.e. the chinoline as well as the backpocket pyridine binding moiety are spread all over 3D space. This can be improved by a modified query that now also contains an additional aromatic feature next to the hinge binder.\n",
    "\n",
    "![4point ph4](images/1py5Ph44pointph4InPocket.png)\n",
    "\n",
    "Running the same alignment as above ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "moePh44Point = \"\"\"#moe:ph4que 2014.09\n",
    "#pharmacophore 7 tag t value *\n",
    "scheme t Unified matchsize i 0 use_Hs i 1 abspos i 0 title t $ useRval i 0 comment s $\n",
    "#feature 4 expr tt color ix x r y r z r r r ebits ix gbits ix m ix\n",
    "Acc df2f2 3.877 7.014 1.448 0.3 0 400 a64cff \n",
    "Acc df2f2 7.22 11.077 5.625 0.3 0 400 a64cff \n",
    "Don f20df2 4.778 8.432 7.805 0.3 0 400 a64cff \n",
    "Aro ff8000 1.56433333333334 7.06399999999999 3.135 1.5 0 400 a64cff\n",
    "#endpharmacophore\"\"\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[12.247548982275546,\n",
       " 11.111090082701395,\n",
       " 3.826678589119382,\n",
       " 11.66292566298035,\n",
       " 8.204876289025478,\n",
       " 0.8101042567665502,\n",
       " 8.125825259923843,\n",
       " 5.6827776196306345,\n",
       " 6.945428752476992,\n",
       " 0.7102500150487856]"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ph4Info = parsePh4Content(moePh44Point)\n",
    "feats,radii = convertFeatureInfo(ph4Info)\n",
    "pcophore = Pharmacophore.Pharmacophore(feats)\n",
    "applyRadiiToBounds(radii,pcophore)\n",
    "ligand = Chem.MolFromSmiles(\"c1ccc(-c2n[nH]cc2-c2ccnc3ccccc23)nc1\")\n",
    "boundsMat = rdDistGeom.GetMoleculeBoundsMatrix(ligand)\n",
    "canMatch,allMatches = EmbedLib.MatchPharmacophoreToMol(ligand,featFactory,pcophore)\n",
    "if canMatch:\n",
    "  failed,boundsMatMatched,matched,matchDetails = EmbedLib.MatchPharmacophore(allMatches,boundsMat,\n",
    "                                                                             pcophore,useDownsampling=True)\n",
    "atomMatch = [list(x.GetAtomIds()) for x in matched]\n",
    "ligH = Chem.AddHs(ligand)\n",
    "bm,embeddings,numFail = EmbedLib.EmbedPharmacophore(ligH,atomMatch,pcophore,count=10)\n",
    "TransformEmbeddings(pcophore,embeddings,atomMatch)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adding constraints into the embedding clearly worsens the SSD values for many alignments. Put differently, even though the distance matrix allows the ph4 to be fit when the actual matching in 3D space happens only few of the matches can be properly embedded onto the pharmacophore. Still, for example alignment 6 shows a low SSD value and the actual alignment looks good.\n",
    "\n",
    "![Alignment 3](images/4PointAlignment6.png)\n",
    "\n",
    "Note that in this match the H-Bond donor is matched to the pyrazol nitrogen without the hydrogen atom. This is due to the definitions in the BaseFeatures.fdef file and is due to the potential tautomerism of the pyrazole group. \n",
    "This is what we would get if we \"cheat\" the second donor atom away (note that the first full match of the pharmacophore feature mappings is used in the embedding). We could also achieve the same by changing the definitions of the BaseFeatures file but this is easier for now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 Acceptor SingleAtomAcceptor (5,)\n",
      "0 Acceptor SingleAtomAcceptor (12,)\n",
      "0 Acceptor SingleAtomAcceptor (19,)\n",
      "1 Acceptor SingleAtomAcceptor (5,)\n",
      "1 Acceptor SingleAtomAcceptor (12,)\n",
      "1 Acceptor SingleAtomAcceptor (19,)\n",
      "2 Donor SingleAtomDonor (5,)\n",
      "2 Donor SingleAtomDonor (6,)\n",
      "3 Aromatic Arom5 (4, 5, 6, 7, 8)\n",
      "3 Aromatic Arom6 (0, 1, 2, 3, 19, 20)\n",
      "3 Aromatic Arom6 (9, 10, 11, 12, 13, 18)\n",
      "3 Aromatic Arom6 (13, 14, 15, 16, 17, 18)\n"
     ]
    }
   ],
   "source": [
    "for (i,match) in enumerate(allMatches):\n",
    "  for f in match:\n",
    "    print(\"%d %s %s %s\"%(i, f.GetFamily(), f.GetType(), f.GetAtomIds()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "allMatches[2] = (allMatches[2][1],)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 Acceptor SingleAtomAcceptor (5,)\n",
      "0 Acceptor SingleAtomAcceptor (12,)\n",
      "0 Acceptor SingleAtomAcceptor (19,)\n",
      "1 Acceptor SingleAtomAcceptor (5,)\n",
      "1 Acceptor SingleAtomAcceptor (12,)\n",
      "1 Acceptor SingleAtomAcceptor (19,)\n",
      "2 Donor SingleAtomDonor (6,)\n",
      "3 Aromatic Arom5 (4, 5, 6, 7, 8)\n",
      "3 Aromatic Arom6 (0, 1, 2, 3, 19, 20)\n",
      "3 Aromatic Arom6 (9, 10, 11, 12, 13, 18)\n",
      "3 Aromatic Arom6 (13, 14, 15, 16, 17, 18)\n"
     ]
    }
   ],
   "source": [
    "for (i,match) in enumerate(allMatches):\n",
    "  for f in match:\n",
    "    print(\"%d %s %s %s\"%(i, f.GetFamily(), f.GetType(), f.GetAtomIds()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.8460812176612649,\n",
       " 11.59755647293786,\n",
       " 2.4969798350293004,\n",
       " 12.165228799867712,\n",
       " 8.098470400659394,\n",
       " 0.701126156123749,\n",
       " 8.488922578961997,\n",
       " 8.821280152742723,\n",
       " 9.94333946662374,\n",
       " 0.637983263851595]"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ligand = Chem.MolFromSmiles(\"c1ccc(-c2n[nH]cc2-c2ccnc3ccccc23)nc1\")\n",
    "boundsMat = rdDistGeom.GetMoleculeBoundsMatrix(ligand)\n",
    "failed,boundsMatMatched,matched,matchDetails = EmbedLib.MatchPharmacophore(allMatches,boundsMat,\n",
    "                                                                           pcophore,useDownsampling=True)\n",
    "atomMatch = [list(x.GetAtomIds()) for x in matched]\n",
    "ligH = Chem.AddHs(ligand)\n",
    "bm,embeddings,numFail = EmbedLib.EmbedPharmacophore(ligH,atomMatch,pcophore,count=10)\n",
    "TransformEmbeddings(pcophore,embeddings,atomMatch)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, alignment 6 shows a good SSD value which is also reflected by a good alignment (this time with the correct Donor alignment).\n",
    "\n",
    "![Donor alignment](images/4PointAlignmentCheat6N.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Caveat !!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the next steps could be to completely \"pin down\" the ph4 so that the compound really fits the ph4. In order to do that, the backpocket pyridine would need to be brought into plane with the rest of the molecule. To achieve that, the aromatic ph4 feature needs to be included. \n",
    "This is where things can get tricky ... in this case the acceptor feature of the pyridine ring system is also part of the aromatic system which means that during the bounds smoothing process an interference between the two will hinder a positive matching. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "moePh45Point = \"\"\"#moe:ph4que 2014.09\n",
    "#pharmacophore 7 tag t value *\n",
    "scheme t Unified matchsize i 0 use_Hs i 1 abspos i 0 title t $ useRval i\n",
    "0 comment s $\n",
    "#feature 5 expr tt color ix x r y r z r r r ebits ix gbits ix m ix\n",
    "Acc df2f2 3.877 7.014 1.448 1 0 400 a64cff \n",
    "Acc df2f2 7.22 11.077 5.625 1 0 400 a64cff \n",
    "Don f20df2 4.778 8.432 7.805 1 0 400 a64cff \n",
    "Aro ff8000 1.56433333333334 7.06399999999999 3.135 1 0 400 a64cff \n",
    "Aro ff8000 6.68983333333333 11.6213333333333 4.498 1 0 400 a64cff\n",
    "#endpharmacophore\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ph4Info = parsePh4Content(moePh45Point)\n",
    "feats,radii = convertFeatureInfo(ph4Info)\n",
    "pcophore = Pharmacophore.Pharmacophore(feats)\n",
    "applyRadiiToBounds(radii,pcophore)\n",
    "ligand = Chem.MolFromSmiles(\"c1ccc(-c2n[nH]cc2-c2ccnc3ccccc23)nc1\")\n",
    "boundsMat = rdDistGeom.GetMoleculeBoundsMatrix(ligand)\n",
    "canMatch,allMatches = EmbedLib.MatchPharmacophoreToMol(ligand,featFactory,pcophore)\n",
    "if canMatch:\n",
    "  failed,boundsMatMatched,matched,matchDetails = EmbedLib.MatchPharmacophore(allMatches,boundsMat,\n",
    "                                                                             pcophore,useDownsampling=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "failed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__When setting up ph4 queries, make sure that multi-atom pharmacophoric features do not include other pharmacophoric features!__ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Running over multiple molecules"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have know how to align the co-crystallised ligand back to its pharmacophore query we can run a couple of molecules that should fit into the pocket of that target enzyme (here Ren et al. doi:10.1016/j.ejmech.2009.07.008). We use the first 5 molecules from the list of actives and use the 4-point pharmacophore used above but with a radius of 1.0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "moePh44Point = \"\"\"#moe:ph4que 2014.09\n",
    "#pharmacophore 7 tag t value *\n",
    "scheme t Unified matchsize i 0 use_Hs i 1 abspos i 0 title t $ useRval i 0 comment s $\n",
    "#feature 4 expr tt color ix x r y r z r r r ebits ix gbits ix m ix\n",
    "Acc df2f2 3.877 7.014 1.448 1.0 0 400 a64cff \n",
    "Acc df2f2 7.22 11.077 5.625 1.0 0 400 a64cff \n",
    "Don f20df2 4.778 8.432 7.805 1.0 0 400 a64cff \n",
    "Aro ff8000 1.56433333333334 7.06399999999999 3.135 1.0 0 400 a64cff\n",
    "#endpharmacophore\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "Couldn't embed molecule 3\n",
      "4\n"
     ]
    }
   ],
   "source": [
    "from operator import itemgetter\n",
    "ph4Info = parsePh4Content(moePh44Point)\n",
    "feats,radii = convertFeatureInfo(ph4Info)\n",
    "pcophore = Pharmacophore.Pharmacophore(feats)\n",
    "applyRadiiToBounds(radii,pcophore)\n",
    "molSmiles = ['Cc1cccc(c2n[nH]cc2c3ccc4ncccc4n3)n1','Cc1cccnc1c2nc(N)sc2c3nc4cccnc4cc3',\n",
    "             'Cc1cccc(c2[nH]c(CNc5cc(C(=O)N)ccc5)nc2c3ccc4nccnc4c3)n1','Clc1cccc(c2nc(N)sc2c3ccc4ncccc4n3)c1',\n",
    "                 'n1ccccc1c2nn3CCCc3c2c4ccnc5cc(NC(=O)NCCN(C)C)ccc45']\n",
    "mols = [Chem.MolFromSmiles(smi) for smi in molSmiles]\n",
    "res = []\n",
    "for i,mol in enumerate(mols):\n",
    "  print(i)\n",
    "  boundsMat = rdDistGeom.GetMoleculeBoundsMatrix(mol)\n",
    "  canMatch,allMatches = EmbedLib.MatchPharmacophoreToMol(mol,featFactory,pcophore)\n",
    "  if canMatch:\n",
    "    failed,boundsMatMatched,matched,matchDetails = EmbedLib.MatchPharmacophore(allMatches,boundsMat,\n",
    "                                                                               pcophore,useDownsampling=True)\n",
    "    if failed:\n",
    "      print(\"Couldn't embed molecule %d\"%i)\n",
    "      continue\n",
    "  else:\n",
    "    print(\"Couldn't match molecule %d\"%i)\n",
    "    continue\n",
    "  atomMatch = [list(x.GetAtomIds()) for x in matched]\n",
    "  try:\n",
    "    molH = Chem.AddHs(mol)\n",
    "    bm,embeddings,numFail = EmbedLib.EmbedPharmacophore(molH,atomMatch,pcophore,count=10)\n",
    "  except ValueError:\n",
    "    print (\"Bounds smoothing failed for molecule %d\"%i)\n",
    "    continue\n",
    "  SSDs = TransformEmbeddings(pcophore,embeddings,atomMatch) \n",
    "  bestFitIndex = min(enumerate(SSDs), key=itemgetter(1))[0] \n",
    "  res.append((SSDs[bestFitIndex],embeddings[bestFitIndex]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(5.4789081986121175, <rdkit.Chem.rdchem.Mol at 0x7f8ed5e93730>),\n",
       " (6.390745025093736, <rdkit.Chem.rdchem.Mol at 0x7f8ed5e930a8>),\n",
       " (5.4231998457926665, <rdkit.Chem.rdchem.Mol at 0x7f8ed5e7ddb8>),\n",
       " (6.303848386994872, <rdkit.Chem.rdchem.Mol at 0x7f8ed5e93368>)]"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As one can see based on RMSD values, the overall alignment to the pharmacophore is not very good. This is also reflected when looking at how the molecules are aligned:\n",
    "\n",
    "![multiMol alignment](images/4PointAlignmentAllN.png)\n",
    "\n",
    "Of course - this doesn't mean that the method itself is not good - it just shows that the molecules are not fitting the pharmacophore. The above example is merely there to show how to run the method on multiples molecules.\n",
    "In addition it shows that one can use the basic alignment values in order to quickly remove molecules that do not properly fit the pharmacophore. However, be cautious since the conformations of the aligned molecules are not really good - but that's for another tutorial ... "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}