       1 #!/bin/env python
       2 # File: RDKitUtil.py
       3 # Author: Manish Sud <msud@san.rr.com>
       4 #
       5 # Copyright (C) 2018 Manish Sud. All rights reserved.
       6 #
       7 # The functionality available in this file is implemented using RDKit, an
       8 # open source toolkit for cheminformatics developed by Greg Landrum.
       9 #
      10 # This file is part of MayaChemTools.
      11 #
      12 # MayaChemTools is free software; you can redistribute it and/or modify it under
      13 # the terms of the GNU Lesser General Public License as published by the Free
      14 # Software Foundation; either version 3 of the License, or (at your option) any
      15 # later version.
      16 #
      17 # MayaChemTools is distributed in the hope that it will be useful, but without
      18 # any warranty; without even the implied warranty of merchantability of fitness
      19 # for a particular purpose.  See the GNU Lesser General Public License for more
      20 # details.
      21 #
      22 # You should have received a copy of the GNU Lesser General Public License
      23 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or
      24 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330,
      25 # Boston, MA, 02111-1307, USA.
      26 #
      27
      28 from __future__ import print_function
      29
      30 import os
      31 import sys
      32 import re
      33
      34 from rdkit import Chem
      35 from rdkit.Chem import AllChem
      36 from rdkit.Chem import Draw
      37
      38 import MiscUtil
      39
      40 __all__ = ["GetInlineSVGForMolecule", "GetInlineSVGForMolecules", "GetMolName", "GetSVGForMolecule", "GetSVGForMolecules", "IsMolEmpty", "ReadMolecules", "ReadAndValidateMolecules", "ReadMoleculesFromSDFile", "ReadMoleculesFromMolFile", "ReadMoleculesFromMol2File", "ReadMoleculesFromPDBFile", "ReadMoleculesFromSMILESFile", "WriteMolecules"]
      41
      42 def GetMolName(Mol, MolNum = None):
      43     """Get molecule name.
      44
      45     Arguments:
      46         Mol (object): RDKit molecule object.
      47         MolNum (int or None): Molecule number in input file.
      48
      49     Returns:
      50         str : Molname corresponding to _Name property of a molecule, generated
      51             from specieid MolNum using the format "Mol%d" % MolNum, or an
      52             empty string.
      53
      54     """
      55
      56     MolName = ''
      57     if Mol.HasProp("_Name"):
      58         MolName = Mol.GetProp("_Name")
      59
      60     if not len(MolName):
      61         if MolNum is not None:
      62             MolName = "Mol%d" % MolNum
      63
      64     return MolName
      65
      66 def GetInlineSVGForMolecule(Mol, Width, Height, Legend = None, AtomListToHighlight = None, BondListToHighlight = None, BoldText = True):
      67     """Get SVG image text for a molecule suitable for inline embedding into a HTML page.
      68
      69     Arguments:
      70         Mol (object): RDKit molecule object.
      71         Width (int): Width of a molecule image in pixels.
      72         Height (int): Height of a molecule image in pixels.
      73         Legend (str): Text to display under the image.
      74         AtomListToHighlight (list): List of atoms to highlight.
      75         BondListToHighlight (list): List of bonds to highlight.
      76         BoldText (bool): Flag to make text bold in the image of molecule.
      77
      78     Returns:
      79         str : SVG image text for inline embedding into a HTML page using "img"
      80             tag: <img src="data:image/svg+xml;charset=UTF-8,SVGImageText>
      81
      82     """
      83
      84     SVGText = GetSVGForMolecule(Mol, Width, Height, Legend, AtomListToHighlight, BondListToHighlight, BoldText)
      85     return _ModifySVGForInlineEmbedding(SVGText)
      86
      87 def GetInlineSVGForMolecules(Mols, MolsPerRow, MolWidth, MolHeight, Legends = None, AtomListsToHighlight = None, BondListsToHighLight = None, BoldText = True):
      88     """Get SVG image text for  molecules suitable for inline embedding into a HTML page.
      89
      90     Arguments:
      91         Mols (list): List of RDKit molecule objects.
      92         MolsPerRow (int): Number of molecules per row.
      93         Width (int): Width of a molecule image in pixels.
      94         Height (int): Height of a molecule image in pixels.
      95         Legends (list): List containing strings to display under images.
      96         AtomListsToHighlight (list): List of lists containing atoms to highlight
      97             for molecules.
      98         BondListsToHighlight (list): List of lists containing bonds to highlight
      99             for molecules
     100         BoldText (bool): Flag to make text bold in the image of molecules.
     101
     102     Returns:
     103         str : SVG image text for inline embedding into a HTML page using "img"
     104             tag: <img src="data:image/svg+xml;charset=UTF-8,SVGImageText>
     105
     106     """
     107
     108     SVGText = GetSVGForMolecules(Mols, MolsPerRow, MolWidth, MolHeight, Legends, AtomListsToHighlight, BondListsToHighLight, BoldText)
     109     return _ModifySVGForInlineEmbedding(SVGText)
     110
     111 def _ModifySVGForInlineEmbedding(SVGText):
     112     """Modify SVG for inline embedding into a HTML page using "img" tag."""
     113
     114     # Take out all tags till the start of '<svg' tag...
     115     Pattern = re.compile("^.*<svg", re.I | re.S)
     116     SVGText = Pattern.sub("<svg", SVGText)
     117
     118     # Add an extra space before the "width=..." tag. Otherwise, inline embedding may
     119     # cause the following XML error on some browsers due to start of the "width=..."
     120     # at the begining of the line in <svg ...> tag:
     121     #
     122     #  XML5607: Whitespace expected.
     123     #
     124     SVGText = re.sub("width='", " width='", SVGText, flags = re.I)
     125
     126     # Take out trailing new line...
     127     SVGText = SVGText.strip()
     128
     129     return SVGText
     130
     131 def GetSVGForMolecule(Mol, Width, Height, Legend = None, AtomListToHighlight = None, BondListToHighlight = None, BoldText = True):
     132     """Get SVG image text for a molecule suitable for viewing in a browser.
     133
     134     Arguments:
     135         Mol (object): RDKit molecule object.
     136         Width (int): Width of a molecule image in pixels.
     137         Height (int): Height of a molecule image in pixels.
     138         Legend (str): Text to display under the image.
     139         AtomListToHighlight (list): List of atoms to highlight.
     140         BondListToHighlight (list): List of bonds to highlight.
     141         BoldText (bool): Flag to make text bold in the image of molecule.
     142
     143     Returns:
     144         str : SVG image text for writing to a SVG file for viewing in a browser.
     145
     146     """
     147
     148     Mols = [Mol]
     149
     150     MolsPerRow = 1
     151     MolWidth = Width
     152     MolHeight = Height
     153
     154     Legends = [Legend] if Legend is not None else None
     155     AtomListsToHighlight = [AtomListToHighlight] if AtomListToHighlight is not None else None
     156     BondListsToHighLight = [BondListsToHighLight] if BondListToHighlight is not None else None
     157
     158     return GetSVGForMolecules(Mols, MolsPerRow, MolWidth, MolHeight, Legends, AtomListsToHighlight, BondListsToHighLight, BoldText)
     159
     160 def GetSVGForMolecules(Mols, MolsPerRow, MolWidth, MolHeight, Legends = None, AtomListsToHighlight = None, BondListsToHighlight = None, BoldText = True):
     161     """Get SVG image text for molecules suitable for viewing in a browser.
     162
     163     Arguments:
     164         Mols (list): List of RDKit molecule objects.
     165         MolsPerRow (int): Number of molecules per row.
     166         Width (int): Width of a molecule image in pixels.
     167         Height (int): Height of a molecule image in pixels.
     168         Legends (list): List containing strings to display under images.
     169         AtomListsToHighlight (list): List of lists containing atoms to highlight
     170             for molecules.
     171         BondListsToHighlight (list): List of lists containing bonds to highlight
     172             for molecules
     173         BoldText (bool): Flag to make text bold in the image of molecules.
     174
     175     Returns:
     176         str : SVG image text for writing to a SVG file for viewing in a browser.
     177
     178     """
     179
     180     SVGText = Draw.MolsToGridImage(Mols, molsPerRow = MolsPerRow, subImgSize = (MolWidth,MolHeight), legends = Legends, highlightAtomLists = AtomListsToHighlight, highlightBondLists = BondListsToHighlight, useSVG = True)
     181
     182     return _ModifySVGForBrowserViewing(SVGText, BoldText)
     183
     184 def _ModifySVGForBrowserViewing(SVGText, BoldText = True):
     185     """Modify SVG for loading into a browser."""
     186
     187     # It appears that the string 'xmlns:svg' needs to be replaced with 'xmlns' in the
     188     # SVG image string generated by RDKit. Otherwise, the image doesn't load in
     189     # web browsers.
     190     #
     191     if re.search("xmlns:svg", SVGText, re.I):
     192         SVGText = re.sub("xmlns:svg", "xmlns", SVGText, flags = re.I)
     193
     194     # Make text bold...
     195     if BoldText:
     196         SVGText = re.sub("font-weight:normal;", "font-weight:bold;", SVGText, flags = re.I)
     197
     198     return SVGText
     199
     200
     201 def IsMolEmpty(Mol):
     202     """Check for the presence of atoms in a molecule.
     203
     204     Arguments:
     205         Mol (object): RDKit molecule object.
     206
     207     Returns:
     208         bool : True - No atoms in molecule; Otherwise, false.
     209
     210     """
     211
     212     Status = False if Mol.GetNumAtoms() else True
     213
     214     return Status
     215
     216 def ReadAndValidateMolecules(FileName, **KeyWordArgs):
     217     """Read molecules from an input file, validate all molecule objects, and return
     218     a list of valid and non-valid molecule objects along with their counts.
     219
     220     Arguments:
     221         FileName (str): Name of a file with complete path.
     222         **KeyWordArgs (dictionary) : Parameter name and value pairs for reading and
     223             processing molecules.
     224
     225     Returns:
     226         list : List of valid RDKit molecule objects.
     227         int : Number of total molecules in input file.
     228         int : Number of valid molecules in input file.
     229
     230     Notes:
     231         The file extension is used to determine type of the file and set up an appropriate
     232         file reader.
     233
     234     """
     235
     236     AllowEmptyMols = True
     237     if "AllowEmptyMols" in KeyWordArgs:
     238         AllowEmptyMols = KeyWordArgs["AllowEmptyMols"]
     239
     240     Mols = ReadMolecules(FileName, **KeyWordArgs)
     241
     242     if AllowEmptyMols:
     243         ValidMols = [Mol for Mol in Mols if Mol is not None]
     244     else:
     245         ValidMols = []
     246         MolCount = 0
     247         for Mol in Mols:
     248             MolCount += 1
     249             if Mol is None:
     250                 continue
     251
     252             if IsMolEmpty(Mol):
     253                 MolName = GetMolName(Mol, MolCount)
     254                 MiscUtil.PrintWarning("Ignoring empty molecule: %s" % MolName)
     255                 continue
     256
     257             ValidMols.append(Mol)
     258
     259     MolCount = len(Mols)
     260     ValidMolCount = len(ValidMols)
     261
     262     return (ValidMols, MolCount, ValidMolCount)
     263
     264 def ReadMolecules(FileName, **KeyWordArgs):
     265     """Read molecules from an input file without performing any validation
     266     and creation of molecule objects.
     267
     268     Arguments:
     269         FileName (str): Name of a file with complete path.
     270         **KeyWordArgs (dictionary) : Parameter name and value pairs for reading and
     271             processing molecules.
     272
     273     Returns:
     274         list : List of RDKit molecule objects.
     275
     276     Notes:
     277         The file extension is used to determine type of the file and set up an appropriate
     278         file reader.
     279
     280     """
     281
     282     # Set default values for possible arguments...
     283     ReaderArgs = {"Sanitize": True, "RemoveHydrogens": True, "StrictParsing": True,  "SMILESDelimiter" : ' ', "SMILESColumn": 1, "SMILESNameColumn": 2, "SMILESTitleLine": True }
     284
     285     # Set specified values for possible arguments...
     286     for Arg in ReaderArgs:
     287         if Arg in KeyWordArgs:
     288             ReaderArgs[Arg] = KeyWordArgs[Arg]
     289
     290     # Modify specific valeus for SMILES...
     291     if MiscUtil.CheckFileExt(FileName, "smi csv tsv txt"):
     292         Args = ["Sanitize", "SMILESTitleLine"]
     293         for Arg in Args:
     294             if ReaderArgs[Arg] is True:
     295                 ReaderArgs[Arg] = 1
     296             else:
     297                 ReaderArgs[Arg] = 0
     298
     299     Mols = []
     300     if MiscUtil.CheckFileExt(FileName, "sdf sd"):
     301         return ReadMoleculesFromSDFile(FileName, ReaderArgs["Sanitize"], ReaderArgs["RemoveHydrogens"], ReaderArgs['StrictParsing'])
     302     elif MiscUtil.CheckFileExt(FileName, "mol"):
     303         return ReadMoleculesFromMolFile(FileName, ReaderArgs["Sanitize"], ReaderArgs["RemoveHydrogens"], ReaderArgs['StrictParsing'])
     304     elif MiscUtil.CheckFileExt(FileName, "mol2"):
     305         return ReadMoleculesFromMol2File(FileName, ReaderArgs["Sanitize"], ReaderArgs["RemoveHydrogens"])
     306     elif MiscUtil.CheckFileExt(FileName, "pdb"):
     307         return ReadMoleculesFromPDBFile(FileName, ReaderArgs["Sanitize"], ReaderArgs["RemoveHydrogens"])
     308     elif MiscUtil.CheckFileExt(FileName, "smi txt csv tsv"):
     309         SMILESColumnIndex = ReaderArgs["SMILESColumn"] - 1
     310         SMILESNameColumnIndex = ReaderArgs["SMILESNameColumn"] - 1
     311         return ReadMoleculesFromSMILESFile(FileName, ReaderArgs["SMILESDelimiter"], SMILESColumnIndex, SMILESNameColumnIndex, ReaderArgs["SMILESTitleLine"], ReaderArgs["Sanitize"])
     312     else:
     313         MiscUtil.PrintWarning("RDKitUtil.ReadMolecules: Non supported file type: %s" % FileName)
     314
     315     return Mols
     316
     317 def ReadMoleculesFromSDFile(FileName, Sanitize = True, RemoveHydrogens = True, StrictParsing = True):
     318     """Read molecules from a SD file.
     319
     320     Arguments:
     321         FileName (str): Name of a file with complete path.
     322         Sanitize (bool): Sanitize molecules.
     323         RemoveHydrogens (bool): Remove hydrogens from molecules.
     324         StrictParsing (bool): Perform strict parsing.
     325
     326     Returns:
     327         list : List of RDKit molecule objects.
     328
     329     """
     330     return  Chem.SDMolSupplier(FileName, sanitize = Sanitize, removeHs = RemoveHydrogens, strictParsing = StrictParsing)
     331
     332 def ReadMoleculesFromMolFile(FileName, Sanitize = True, RemoveHydrogens = True, StrictParsing = True):
     333     """Read molecule from a MDL Mol file.
     334
     335     Arguments:
     336         FileName (str): Name of a file with complete path.
     337         Sanitize (bool): Sanitize molecules.
     338         RemoveHydrogens (bool): Remove hydrogens from molecules.
     339         StrictParsing (bool): Perform strict parsing.
     340
     341     Returns:
     342         list : List of RDKit molecule objects.
     343
     344     """
     345
     346     Mols = []
     347     Mols.append(Chem.MolFromMolFile(FileName, sanitize = Sanitize, removeHs = RemoveHydrogens, strictParsing = StrictParsing))
     348     return Mols
     349
     350 def ReadMoleculesFromMol2File(FileName, Sanitize = True, RemoveHydrogens = True):
     351     """Read molecule from a Tripos Mol2  file.
     352
     353     Arguments:
     354         FileName (str): Name of a file with complete path.
     355         Sanitize (bool): Sanitize molecules.
     356         RemoveHydrogens (bool): Remove hydrogens from molecules.
     357
     358     Returns:
     359         list : List of RDKit molecule objects.
     360
     361     """
     362
     363     Mols = []
     364     Mols.append(Chem.MolFromMol2File(FileName,  sanitize = Sanitize, removeHs = RemoveHydrogens))
     365     return Mols
     366
     367 def ReadMoleculesFromPDBFile(FileName, Sanitize = True, RemoveHydrogens = True):
     368     """Read molecule from a PDB  file.
     369
     370     Arguments:
     371         FileName (str): Name of a file with complete path.
     372         Sanitize (bool): Sanitize molecules.
     373         RemoveHydrogens (bool): Remove hydrogens from molecules.
     374
     375     Returns:
     376         list : List of RDKit molecule objects.
     377
     378     """
     379
     380     Mols = []
     381     Mols.append(Chem.MolFromPDBFile(FileName,  sanitize = Sanitize, removeHs = RemoveHydrogens))
     382     return Mols
     383
     384 def ReadMoleculesFromSMILESFile(FileName, SMILESDelimiter = ' ', SMILESColIndex = 0, SMILESNameColIndex = 1, SMILESTitleLine = 1, Sanitize = 1):
     385     """Read molecules from a SMILES file.
     386
     387     Arguments:
     388         SMILESDelimiter (str): Delimiter for parsing SMILES line
     389         SMILESColIndex (int): Column index containing SMILES string.
     390         SMILESNameColIndex (int): Column index containing molecule name.
     391         SMILESTitleLine (int): Flag to indicate presence of title line.
     392         Sanitize (int): Sanitize molecules.
     393
     394     Returns:
     395         list : List of RDKit molecule objects.
     396
     397     """
     398
     399     return  Chem.SmilesMolSupplier(FileName, delimiter = SMILESDelimiter, smilesColumn = SMILESColIndex, nameColumn = SMILESNameColIndex, titleLine = SMILESTitleLine, sanitize = Sanitize)
     400
     401 def MoleculesWriter(FileName, **KeyWordArgs):
     402     """Set up a molecule writer.
     403
     404     Arguments:
     405         FileName (str): Name of a file with complete path.
     406         **KeyWordArgs (dictionary) : Parameter name and value pairs for writing and
     407             processing molecules.
     408
     409     Returns:
     410         RDKit object : Molecule writer.
     411
     412     Notes:
     413         The file extension is used to determine type of the file and set up an appropriate
     414         file writer.
     415
     416     """
     417
     418     # Set default values for possible arguments...
     419     WriterArgs = {"Compute2DCoords" : False, "Kekulize": False, "SMILESDelimiter" : ' ', "SMILESIsomeric": True, "SMILESTitleLine": True }
     420
     421     # Set specified values for possible arguments...
     422     for Arg in WriterArgs:
     423         if Arg in KeyWordArgs:
     424             WriterArgs[Arg] = KeyWordArgs[Arg]
     425
     426     Writer = None
     427     if MiscUtil.CheckFileExt(FileName, "sdf sd"):
     428         Writer = Chem.SDWriter(FileName)
     429         if WriterArgs["Kekulize"]:
     430             Writer.SetKekulize(True)
     431     elif MiscUtil.CheckFileExt(FileName, "pdb"):
     432         Writer = Chem.PDBWriter(FileName)
     433     elif MiscUtil.CheckFileExt(FileName, "smi"):
     434         # Text for the name column in the title line. Blank indicates not to include name column
     435         # in the output file...
     436         NameHeader = 'Name'
     437         Writer = Chem.SmilesWriter(FileName, delimiter = WriterArgs["SMILESDelimiter"], nameHeader = NameHeader, includeHeader = WriterArgs["SMILESTitleLine"],  isomericSmiles = WriterArgs["SMILESIsomeric"], kekuleSmiles = WriterArgs["Kekulize"])
     438     else:
     439         MiscUtil.PrintWarning("RDKitUtil.WriteMolecules: Non supported file type: %s" % FileName)
     440
     441     return Writer
     442
     443 def WriteMolecules(FileName, Mols, **KeyWordArgs):
     444     """Write molecules to an output file.
     445
     446     Arguments:
     447         FileName (str): Name of a file with complete path.
     448         Mols (list): List of RDKit molecule objects.
     449         **KeyWordArgs (dictionary) : Parameter name and value pairs for writing and
     450             processing molecules.
     451
     452     Returns:
     453         int : Number of total molecules.
     454         int : Number of processed molecules written to output file.
     455
     456     Notes:
     457         The file extension is used to determine type of the file and set up an appropriate
     458         file writer.
     459
     460     """
     461
     462     Compute2DCoords = False
     463     if "Compute2DCoords" in KeyWordArgs:
     464         Compute2DCoords = KeyWordArgs["Compute2DCoords"]
     465
     466     MolCount = len(Mols)
     467     ProcessedMolCount = 0
     468
     469     Writer = MoleculesWriter(FileName, **KeyWordArgs)
     470
     471     if Writer is None:
     472         return (MolCount, ProcessedMolCount)
     473
     474     for Mol in Mols:
     475         if Mol is None:
     476             continue
     477
     478         ProcessedMolCount += 1
     479         if Compute2DCoords:
     480             AllChem.Compute2DCoords(Mol)
     481
     482         Writer.write(Mol)
     483
     484     Writer.close()
     485
     486     return (MolCount, ProcessedMolCount)