{"cells":[{"metadata":{},"cell_type":"markdown","source":"## Purpose\n\nA quick notebook to show how to extract key structural information from molecule objects in RDKIt."},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Installing required packages. This can take a while ...')\n# We use a specific version of RDKit to ensure consistent results.\n!conda install -q -y -c rdkit rdkit=2020_03_6\n!pip -q install cairosvg==2.5.2\nprint('DONE.')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Imports"},{"metadata":{"trusted":true},"cell_type":"code","source":"from io import BytesIO\n\nimport numpy as np\nfrom PIL import Image\n\nfrom rdkit import Chem\nfrom rdkit.Chem import AllChem, Draw\nfrom rdkit.Chem.Draw import rdMolDraw2D\n\nimport lxml.etree as et\nimport cairosvg\n\nfrom matplotlib import pyplot as plt","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Implementation\n\nBelow, the atoms and bonds in the molecule are shown with information about each atom:\n\n* the atom's element symbol\n* the index of the atom within the RDKiet molecule object\n* the number of hydrogens connected to the atom\n* the coordinates of the atom in the rendered output image\n* the type of chirality for the atom (if applicable)\n* whether the atom is part of a ring (shown as '⭕')\n\nand its bonds:\n\n* the atoms at each end of the bond (specified by index)\n* bond strength\n* the bond type (this seems only relevant to stereochemisty)\n* whether the bond is part of a ring (shown as '⭕')"},{"metadata":{"trusted":true},"cell_type":"code","source":"def svg_to_image(svg, convert_to_greyscale=True):\n    svg_str = et.tostring(svg)\n    # TODO: would prefer to convert SVG dirrectly to a numpy array.\n    png = cairosvg.svg2png(bytestring=svg_str)\n    image = np.array(Image.open(BytesIO(png)), dtype=np.float32)\n    # Naive greyscale conversion.\n    if convert_to_greyscale:\n        image = image.mean(axis=-1)\n    return image\n\n\ndef drawmol(mol, size, add_hs=False):\n    if add_hs:\n        mol = Chem.AddHs(mol)\n    d = Draw.rdMolDraw2D.MolDraw2DSVG(2*size, size)\n    AllChem.Compute2DCoords(mol)\n    mol = rdMolDraw2D.PrepareMolForDrawing(mol) \n    d.DrawMolecule(mol)\n    d.FinishDrawing()\n    drawing =  d.GetDrawingText()\n    svg = et.fromstring(drawing.encode('iso-8859-1'))\n    img = 1 - svg_to_image(svg)/255\n    height, width = img.shape\n    xs, ys = [], []\n    atoms = mol.GetAtoms()\n    for iatom in range(mol.GetNumAtoms()):\n        p = d.GetDrawCoords(iatom)\n        x, y = p.x, p.y\n        xs.append(x) ; ys.append(y)\n        atom = atoms[iatom]\n        num_h = atom.GetNumExplicitHs() + atom.GetNumImplicitHs()\n        chiral_tag = atom.GetChiralTag()\n        chiral_label = chiral_tag != Chem.rdchem.ChiralType.CHI_UNSPECIFIED and f' {chiral_tag}' or ''\n        ring_atom = atom.IsInRing()\n        print(f'{atom.GetSymbol()} [#{iatom}]: hydrogens={num_h}, coords=({x:0.1f}, {y:0.1f}){chiral_label}{ring_atom and \" ⭕\" or \"\"}')\n        for ibond, bond in enumerate(atom.GetBonds()):\n            stereo = bond.GetStereo()\n            stereo_label = stereo != Chem.rdchem.BondStereo.STEREONONE and f' {stereo}' or ''\n            bond_direction = bond.GetBondDir()\n            bond_direction_label = Chem.rdchem.BondDir.NONE != bond_direction and f' {bond_direction} →' or ' →'\n            in_ring_label = bond.IsInRing() and ' ⭕' or ''\n            bond_type = bond.GetBondType()\n            bond_strength = bond.GetBondTypeAsDouble()\n            from_atom = bond.GetBeginAtomIdx()\n            to_atom = bond.GetEndAtomIdx()\n            print(f'\\t[#{from_atom}] → {bond_type} ({bond_strength}){bond_direction_label} [#{to_atom}]{stereo_label}{in_ring_label}')\n    display(mol)\n    coords = np.array([ys, xs])\n    coords = coords.round().astype(np.int32)\n    coords[0] = np.clip(coords[0], 0, height - 1)\n    coords[1] = np.clip(coords[1], 0, width - 1)\n    atom_layer = np.zeros_like(img)\n    atom_layer[coords[0], coords[1]] = 1\n    def smear_layer(layer):\n        layer = layer + np.roll(layer, 1, axis=0)\n        layer = layer + np.roll(layer, -1, axis=0)\n        layer = layer + np.roll(layer, -1, axis=1)\n        layer = layer + np.roll(layer, 1, axis=1)\n        return layer/layer.max()\n    for _ in range(20):\n        atom_layer = smear_layer(atom_layer)\n    img_rgb = np.stack([atom_layer, 0.8*img, atom_layer], axis=-1)\n    plt.figure(figsize=(14, 14)) ; plt.imshow(img_rgb, cmap='gray_r', interpolation='bilinear') ; plt.show()\n    return 'Done.'\n\n\nINCHI = 'InChI=1S/C19H18ClN3O5S/c1-8-11(12(22-28-8)9-6-4-5-7-10(9)20)15(24)21-13-16(25)23-14(18(26)27)19(2,3)29-17(13)23/h4-7,13-14,17H,1-3H3,(H,21,24)(H,26,27)/t13-,14+,17-/m1/s1'\nmol = Chem.inchi.MolFromInchi(INCHI)\ndrawmol(mol, size=200)\n","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}