       1 #!/bin/env python
       2 #
       3 # File: RDKitGenerateMolecularFrameworks.py
       4 # Author: Manish Sud <msud@san.rr.com>
       5 #
       6 # Copyright (C) 2018 Manish Sud. All rights reserved.
       7 #
       8 # The functionality available in this script is implemented using RDKit, an
       9 # open source toolkit for cheminformatics developed by Greg Landrum.
      10 #
      11 # This file is part of MayaChemTools.
      12 #
      13 # MayaChemTools is free software; you can redistribute it and/or modify it under
      14 # the terms of the GNU Lesser General Public License as published by the Free
      15 # Software Foundation; either version 3 of the License, or (at your option) any
      16 # later version.
      17 #
      18 # MayaChemTools is distributed in the hope that it will be useful, but without
      19 # any warranty; without even the implied warranty of merchantability of fitness
      20 # for a particular purpose.  See the GNU Lesser General Public License for more
      21 # details.
      22 #
      23 # You should have received a copy of the GNU Lesser General Public License
      24 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or
      25 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330,
      26 # Boston, MA, 02111-1307, USA.
      27 #
      28
      29 from __future__ import print_function
      30
      31 # Add local python path to the global path and import standard library modules...
      32 import os
      33 import sys;  sys.path.insert(0, os.path.join(os.path.dirname(sys.argv[0]), "..", "lib", "Python"))
      34 import time
      35 import re
      36 import csv
      37
      38 # RDKit imports...
      39 try:
      40     from rdkit import rdBase
      41     from rdkit import Chem
      42     from rdkit.Chem import AllChem
      43     from rdkit.Chem.Scaffolds import MurckoScaffold
      44 except ImportError as ErrMsg:
      45     sys.stderr.write("\nFailed to import RDKit module/package: %s\n" % ErrMsg)
      46     sys.stderr.write("Check/update your RDKit environment and try again.\n\n")
      47     sys.exit(1)
      48
      49 # MayaChemTools imports...
      50 try:
      51     from docopt import docopt
      52     import MiscUtil
      53     import RDKitUtil
      54 except ImportError as ErrMsg:
      55     sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg)
      56     sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n")
      57     sys.exit(1)
      58
      59 ScriptName = os.path.basename(sys.argv[0])
      60 Options = {}
      61 OptionsInfo = {}
      62
      63 def main():
      64     """Start execution of the script"""
      65
      66     MiscUtil.PrintInfo("\n%s (RDK v%s; %s): Starting...\n" % (ScriptName, rdBase.rdkitVersion, time.asctime()))
      67
      68     (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime()
      69
      70     # Retrieve command line arguments and options...
      71     RetrieveOptions()
      72
      73     # Process and validate command line arguments and options...
      74     ProcessOptions()
      75
      76     # Perform actions required by the script...
      77     GenerateMolecularFrameworks()
      78
      79     MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
      80     MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
      81
      82 def GenerateMolecularFrameworks():
      83     """Generate Bemis Murcko molecular framworks."""
      84
      85     Infile = OptionsInfo["Infile"]
      86     Outfile = OptionsInfo["Outfile"]
      87
      88     UseChirality = OptionsInfo["UseChirality"]
      89
      90     RemoveDuplicateFrameworks = OptionsInfo["RemoveDuplicateFrameworks"]
      91     UseGraphFrameworks = OptionsInfo["UseGraphFrameworks"]
      92
      93     SortFrameworks = OptionsInfo["SortFrameworks"]
      94     if SortFrameworks:
      95         FrameworkMolIDs = []
      96         FrameworkMolIDToMolMap = {}
      97         FrameworkMolIDToAtomCountMap = {}
      98
      99         DuplicateFrameworkMolIDs = []
     100         DuplicateFrameworkMolIDToMolMap = {}
     101         DuplicateFrameworkMolIDToAtomCountMap = {}
     102
     103     DuplicatesOutfile = ""
     104     if RemoveDuplicateFrameworks:
     105         DuplicatesOutfile = OptionsInfo["DuplicatesOutfile"]
     106
     107     # Setup a molecule reader...
     108     MiscUtil.PrintInfo("\nProcessing file %s..." % Infile)
     109     Mols  = RDKitUtil.ReadMolecules(Infile, **OptionsInfo["InfileParams"])
     110
     111     # Set up a molecular framework  writer...
     112     Writer = RDKitUtil.MoleculesWriter(Outfile, **OptionsInfo["OutfileParams"])
     113     if Writer is None:
     114         MiscUtil.PrintError("Failed to setup a writer for output fie %s " % Outfile)
     115
     116     # Set up a duplicate molecular framework writer...
     117     if RemoveDuplicateFrameworks:
     118         DuplicatesWriter = RDKitUtil.MoleculesWriter(DuplicatesOutfile, **OptionsInfo["OutfileParams"])
     119         if Writer is None:
     120             MiscUtil.PrintError("Failed to setup a writer for duplicates output fie %s " % DuplicatesOutfile)
     121
     122     if RemoveDuplicateFrameworks:
     123         MiscUtil.PrintInfo("Generating files: %s and %s..." % (Outfile, DuplicatesOutfile))
     124     else:
     125         MiscUtil.PrintInfo("Generating file %s..." % Outfile)
     126
     127     # Process molecules...
     128     MolCount = 0
     129     ValidMolCount = 0
     130
     131     FrameworksCount = 0
     132     UniqueFrameworksCount = 0
     133     DuplicateFrameworksCount = 0
     134
     135     CanonicalSMILESMap = {}
     136
     137     Compute2DCoords = OptionsInfo["OutfileParams"]["Compute2DCoords"]
     138
     139     for Mol in Mols:
     140         MolCount += 1
     141
     142         if Mol is None:
     143             continue
     144
     145         if RDKitUtil.IsMolEmpty(Mol):
     146             MolName = RDKitUtil.GetMolName(Mol, MolCount)
     147             MiscUtil.PrintWarning("Ignoring empty molecule: %s" % MolName)
     148             continue
     149
     150         ValidMolCount += 1
     151
     152         if UseGraphFrameworks:
     153             FrameworksMol = MurckoScaffold.MakeScaffoldGeneric(Mol)
     154         else:
     155             FrameworksMol = MurckoScaffold.GetScaffoldForMol(Mol)
     156
     157         if Compute2DCoords:
     158             AllChem.Compute2DCoords(FrameworksMol)
     159
     160         if SortFrameworks:
     161             HeavyAtomCount = FrameworksMol.GetNumHeavyAtoms()
     162
     163         FrameworksCount += 1
     164
     165         if RemoveDuplicateFrameworks:
     166             CanonicalSMILES = Chem.MolToSmiles(FrameworksMol, isomericSmiles = UseChirality, canonical = True)
     167             if CanonicalSMILES in CanonicalSMILESMap:
     168                 DuplicateFrameworksCount += 1
     169                 if SortFrameworks:
     170                     # Track duplicate frameworks...
     171                     DuplicateFrameworkMolIDs.append(DuplicateFrameworksCount)
     172                     DuplicateFrameworkMolIDToMolMap[DuplicateFrameworksCount] = FrameworksMol
     173                     DuplicateFrameworkMolIDToAtomCountMap[DuplicateFrameworksCount] = HeavyAtomCount
     174                 else:
     175                     # Write it out...
     176                     DuplicatesWriter.write(FrameworksMol)
     177             else:
     178                 UniqueFrameworksCount += 1
     179                 CanonicalSMILESMap[CanonicalSMILES] = CanonicalSMILES
     180                 if SortFrameworks:
     181                     # Track unique frameworks...
     182                     FrameworkMolIDs.append(UniqueFrameworksCount)
     183                     FrameworkMolIDToMolMap[UniqueFrameworksCount] = FrameworksMol
     184                     FrameworkMolIDToAtomCountMap[UniqueFrameworksCount] = HeavyAtomCount
     185                 else:
     186                     # Write it out...
     187                     Writer.write(FrameworksMol)
     188         elif SortFrameworks:
     189             # Track for sorting...
     190             FrameworkMolIDs.append(FrameworksCount)
     191             FrameworkMolIDToMolMap[FrameworksCount] = FrameworksMol
     192             FrameworkMolIDToAtomCountMap[FrameworksCount] = HeavyAtomCount
     193         else:
     194             # Write it out...
     195             Writer.write(FrameworksMol)
     196
     197     if SortFrameworks:
     198         ReverseOrder = OptionsInfo["DescendingSortOrder"]
     199         SortAndWriteFrameworks(Writer, FrameworkMolIDs, FrameworkMolIDToMolMap, FrameworkMolIDToAtomCountMap, ReverseOrder)
     200         if RemoveDuplicateFrameworks:
     201             SortAndWriteFrameworks(DuplicatesWriter, DuplicateFrameworkMolIDs, DuplicateFrameworkMolIDToMolMap, DuplicateFrameworkMolIDToAtomCountMap, ReverseOrder)
     202
     203     Writer.close()
     204     if RemoveDuplicateFrameworks:
     205         DuplicatesWriter.close()
     206
     207     MiscUtil.PrintInfo("\nTotal number of molecular frameworks: %d" % FrameworksCount)
     208     if RemoveDuplicateFrameworks:
     209         MiscUtil.PrintInfo("Number of unique molecular frameworks: %d" % UniqueFrameworksCount)
     210         MiscUtil.PrintInfo("Number of duplicate molecular frameworks: %d" % DuplicateFrameworksCount)
     211
     212     MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount)
     213     MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount)
     214     MiscUtil.PrintInfo("Number of ignored molecules: %d" % (MolCount - ValidMolCount))
     215
     216 def SortAndWriteFrameworks(MolWriter, MolIDs, MolIDToMolMap, MolIDToAtomCountMap, ReverseOrder):
     217     """Sort frameworks and write them out."""
     218     SortedMolIDs = sorted(MolIDs, key = lambda MolID: MolIDToAtomCountMap[MolID], reverse = ReverseOrder)
     219     for MolID in SortedMolIDs:
     220         FrameworksMol = MolIDToMolMap[MolID]
     221         MolWriter.write(FrameworksMol)
     222
     223 def ProcessOptions():
     224     """Process and validate command line arguments and options"""
     225
     226     MiscUtil.PrintInfo("Processing options...")
     227
     228     # Validate options...
     229     ValidateOptions()
     230
     231     OptionsInfo["Infile"] = Options["--infile"]
     232     OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters("--infileParams", Options["--infileParams"], Options["--infile"])
     233
     234     OptionsInfo["Outfile"] = Options["--outfile"]
     235     OptionsInfo["OutfileParams"] = MiscUtil.ProcessOptionOutfileParameters("--outfileParams", Options["--outfileParams"], Options["--infile"], Options["--outfile"])
     236
     237     OptionsInfo["Overwrite"] = Options["--overwrite"]
     238
     239     OptionsInfo["Mode"] = Options["--mode"]
     240     OptionsInfo["UseGraphFrameworks"] = False
     241     if re.match("^GraphFrameworks$", OptionsInfo["Mode"], re.I):
     242         OptionsInfo["UseGraphFrameworks"] = True
     243
     244     OptionsInfo["RemoveDuplicates"] = Options["--removeDuplicates"]
     245     OptionsInfo["RemoveDuplicateFrameworks"] = False
     246     if re.match("^Yes$", OptionsInfo["RemoveDuplicates"], re.I):
     247         OptionsInfo["RemoveDuplicateFrameworks"] = True
     248
     249     # Setup outfile for writing out duplicates...
     250     OptionsInfo["DuplicatesOutfile"] = ""
     251     if OptionsInfo["RemoveDuplicateFrameworks"]:
     252         FileDir, FileName, FileExt = MiscUtil.ParseFileName(OptionsInfo["Outfile"])
     253         OptionsInfo["DuplicatesOutfile"] = "%sDuplicates.%s" % (FileName, FileExt)
     254
     255     OptionsInfo["Sort"] = Options["--sort"]
     256     OptionsInfo["SortFrameworks"] = False
     257     if re.match("^Yes$", OptionsInfo["Sort"], re.I):
     258         OptionsInfo["SortFrameworks"] = True
     259
     260     OptionsInfo["SortOrder"] = Options["--sortOrder"]
     261     OptionsInfo["DescendingSortOrder"] = False
     262     if re.match("^Descending$", OptionsInfo["SortOrder"], re.I):
     263         OptionsInfo["DescendingSortOrder"] = True
     264
     265     OptionsInfo["UseChirality"] = False
     266     if re.match("^yes$", Options["--useChirality"], re.I):
     267         OptionsInfo["UseChirality"] = True
     268
     269 def RetrieveOptions():
     270     """Retrieve command line arguments and options"""
     271
     272     # Get options...
     273     global Options
     274     Options = docopt(_docoptUsage_)
     275
     276     # Set current working directory to the specified directory...
     277     WorkingDir = Options["--workingdir"]
     278     if WorkingDir:
     279         os.chdir(WorkingDir)
     280
     281     # Handle examples option...
     282     if "--examples" in Options and Options["--examples"]:
     283         MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
     284         sys.exit(0)
     285
     286 def ValidateOptions():
     287     """Validate option values"""
     288
     289     MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
     290     MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "sdf sd smi txt csv tsv")
     291
     292     MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "sdf sd smi")
     293     MiscUtil.ValidateOptionsOutputFileOverwrite("-o, --outfile", Options["--outfile"], "--overwrite", Options["--overwrite"])
     294     MiscUtil.ValidateOptionsDistinctFileNames("-i, --infile", Options["--infile"], "-o, --outfile", Options["--outfile"])
     295
     296     MiscUtil.ValidateOptionTextValue("-m, --mode", Options["--mode"], "GraphFrameworks AtomicFrameworks")
     297
     298     MiscUtil.ValidateOptionTextValue("-r, --removeDuplicates", Options["--removeDuplicates"], "yes no")
     299     MiscUtil.ValidateOptionTextValue("-s, --sort", Options["--sort"], "yes no")
     300     MiscUtil.ValidateOptionTextValue("--sortOrder", Options["--sortOrder"], "ascending descending")
     301
     302     MiscUtil.ValidateOptionTextValue("--useChirality", Options["--useChirality"], "yes no")
     303
     304 # Setup a usage string for docopt...
     305 _docoptUsage_ = """
     306 RDKitGenerateMolecularFrameworks.py - Generate Bemis Murcko molecular frameworks
     307
     308 Usage:
     309     RDKitGenerateMolecularFrameworks.py [--infileParams <Name,Value,...>]
     310                                          [--mode <GraphFrameworks or AtomicFrameworks> ]
     311                                          [ --outfileParams <Name,Value,...> ]  [--overwrite] [--removeDuplicates <yes or no>]
     312                                          [--sort <yes or no>] [--sortOrder <ascending or descending>]
     313                                          [--useChirality <yes or no>] [-w <dir>] -i <infile> -o <outfile>
     314     RDKitGenerateMolecularFrameworks.py -h | --help | -e | --examples
     315
     316 Description:
     317     Generate Bemis Murcko [ Ref 133 ] molecular frameworks for molecules. Two types of molecular
     318     frameworks can be generated: Graph or atomic frameworks. The graph molecular framework
     319     is a generic framework. The atom type, hybridization, and bond order is ignore during its
     320     generation. All atoms are set to carbon atoms and all bonds are single bonds. The atom type,
     321     hybridization, and bond order is preserved during generation of atomic molecular frameworks.
     322
     323     The supported input file formats are: SD (.sdf, .sd), SMILES (.smi, .csv, .tsv, .txt)
     324
     325     The supported output file formats are: SD (.sdf, .sd), SMILES (.smi)
     326
     327 Options:
     328     -e, --examples
     329         Print examples.
     330     -h, --help
     331         Print this help message.
     332     -i, --infile <infile>
     333         Input file name.
     334     --infileParams <Name,Value,...>  [default: auto]
     335         A comma delimited list of parameter name and value pairs for reading
     336         molecules from files. The supported parameter names for different file
     337         formats, along with their default values, are shown below:
     338
     339             SD: removeHydrogens,yes,sanitize,yes,strictParsing,yes
     340             SMILES: smilesColumn,1,smilesNameColumn,2,smilesDelimiter,space,
     341                 smilesTitleLine,auto,sanitize,yes
     342
     343         Possible values for smilesDelimiter: space, comma or tab.
     344     -m, --mode <GraphFrameworks or AtomicFrameworks>  [default: GraphFrameworks]
     345         Type of molecular frameworks to generate for molecules. Possible values: GraphFrameworks
     346          or AtomicFrameworks. The graph molecular framework is a generic framework. The atom type,
     347          hybridization, and bond order is ignore during its generation. All atoms are set to carbon atoms
     348          and all bonds are single bonds. The atom type, hybridization, and bond order is preserved
     349          during the generation of atomic molecular frameworks.
     350     -o, --outfile <outfile>
     351         Output file name.
     352     --outfileParams <Name,Value,...>  [default: auto]
     353         A comma delimited list of parameter name and value pairs for writing
     354         molecules to files. The supported parameter names for different file
     355         formats, along with their default values, are shown below:
     356
     357             SD: compute2DCoords,auto
     358             SMILES: smilesDelimiter,space,smilesTitleLine,yes
     359
     360         Default value for compute2DCoords: yes for SMILES input file; no for all other
     361         file types.
     362     --overwrite
     363         Overwrite existing files.
     364     -r, --removeDuplicates <yes or no>  [default: no]
     365         Remove duplicate molecular frameworks. Possible values: yes or no. The duplicate
     366         molecular franworks are identified using canonical SMILES. The removed frameworks
     367         are written to a separate output file.
     368     -s, --sort <yes or no>  [default: no]
     369         Sort molecular frameworks by heavy atom count. Possible values: yes or no.
     370     --sortOrder <ascending or descending>  [default: ascending]
     371         Sorting order for molecular frameworks. Possible values: ascending or descending.
     372     -u, --useChirality <yes or no>  [default: yes]
     373         Use stereochemistry for generation of canonical SMILES strings to identify
     374         duplicate molecular frameworks.
     375     -w, --workingdir <dir>
     376         Location of working directory which defaults to the current directory.
     377
     378 Examples:
     379     To generate graph molecular framworks for molecules and write out a SMILES file,
     380     type:
     381
     382         % RDKitGenerateMolecularFrameworks.py -i Sample.smi -o SampleOut.smi
     383
     384     To generate graph molecular framworks, remove duplicate frameworks for molecules
     385     and write out SD files for unique and duplicate frameworks, type:
     386
     387         % RDKitGenerateMolecularFrameworks.py -m GraphFrameworks -r yes
     388           -i Sample.sdf -o SampleOut.sdf
     389
     390     To generate atomic molecular framworks, remove duplicate frameworks, sort
     391     framworks by heavy atom count in ascending order, write out SMILES files for
     392     unique and duplicate frameworks, type:
     393
     394         % RDKitGenerateMolecularFrameworks.py -m AtomicFrameworks -r yes
     395           -s yes -i Sample.smi -o SampleOut.smi
     396
     397     To generate graph molecular framworks for molecules in a CSV SMILES file,
     398     SMILES strings in column 1, name in olumn 2, emove duplicate frameworks,
     399     sort framworks by heavy atom count in decending order and write out a SD
     400     file, type:
     401
     402         % RDKitGenerateMolecularFrameworks.py -m AtomicFrameworks
     403           --removeDuplicates yes -s yes --sortOrder descending --infileParams
     404           "smilesDelimiter,comma,smilesTitleLine,yes,smilesColumn,1,
     405           smilesNameColumn,2" --outfileParams "compute2DCoords,yes"
     406           -i SampleSMILES.csv -o SampleOut.sdf
     407
     408 Author:
     409     Manish Sud(msud@san.rr.com)
     410
     411 See also:
     412     RDKitConvertFileFormat.py, RDKitDrawMolecules.py, RDKitSearchFunctionalGroups.py,
     413     RDKitSearchSMARTS.py
     414
     415 Copyright:
     416     Copyright (C) 2018 Manish Sud. All rights reserved.
     417
     418     The functionality available in this script is implemented using RDKit, an
     419     open source toolkit for cheminformatics developed by Greg Landrum.
     420
     421     This file is part of MayaChemTools.
     422
     423     MayaChemTools is free software; you can redistribute it and/or modify it under
     424     the terms of the GNU Lesser General Public License as published by the Free
     425     Software Foundation; either version 3 of the License, or (at your option) any
     426     later version.
     427
     428 """
     429
     430 if __name__ == "__main__":
     431     main()