       1 #!/bin/env python
       2 #
       3 # File: RDKitCalculatePartialCharges.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
      37 # RDKit imports...
      38 try:
      39     from rdkit import rdBase
      40     from rdkit import Chem
      41     from rdkit.Chem import AllChem
      42     from rdkit.Chem import rdPartialCharges
      43 except ImportError as ErrMsg:
      44     sys.stderr.write("\nFailed to import RDKit module/package: %s\n" % ErrMsg)
      45     sys.stderr.write("Check/update your RDKit environment and try again.\n\n")
      46     sys.exit(1)
      47
      48 # RDKit dependency imports...
      49 import numpy
      50
      51 # MayaChemTools imports...
      52 try:
      53     from docopt import docopt
      54     import MiscUtil
      55     import RDKitUtil
      56 except ImportError as ErrMsg:
      57     sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg)
      58     sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n")
      59     sys.exit(1)
      60
      61 ScriptName = os.path.basename(sys.argv[0])
      62 Options = {}
      63 OptionsInfo = {}
      64
      65 DescriptorNamesMap = {}
      66
      67 def main():
      68     """Start execution of the script"""
      69
      70     MiscUtil.PrintInfo("\n%s (RDK v%s; %s): Starting...\n" % (ScriptName, rdBase.rdkitVersion, time.asctime()))
      71
      72     (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime()
      73
      74     # Retrieve command line arguments and options...
      75     RetrieveOptions()
      76
      77     # Process and validate command line arguments and options...
      78     ProcessOptions()
      79
      80     # Perform actions required by the script...
      81     CalculatePartialCharges()
      82
      83     MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
      84     MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
      85
      86 def CalculatePartialCharges():
      87     """Calculate partial atomic charges."""
      88
      89     Infile = OptionsInfo["Infile"]
      90     Outfile = OptionsInfo["Outfile"]
      91
      92     MiscUtil.PrintInfo("Calculating partial atomic charges...")
      93
      94     # Setup a molecule reader...
      95     MiscUtil.PrintInfo("\nProcessing file %s..." % Infile)
      96     Mols = RDKitUtil.ReadMolecules(Infile, **OptionsInfo["InfileParams"])
      97
      98     # Setup a writer...
      99     Compute2DCoords = OptionsInfo["OutfileParams"]["Compute2DCoords"]
     100     Writer = RDKitUtil.MoleculesWriter(Outfile, **OptionsInfo["OutfileParams"])
     101
     102     if Writer is None:
     103         MiscUtil.PrintError("Failed to setup a writer for output fie %s " % Outfile)
     104
     105     MiscUtil.PrintInfo("Generating file %s..." % Outfile)
     106
     107     # Process molecules...
     108     MolCount, ValidMolCount, CalcFailedCount = [0] * 3
     109
     110     for Mol in Mols:
     111         MolCount += 1
     112         if Mol is None:
     113             continue
     114
     115         if RDKitUtil.IsMolEmpty(Mol):
     116             MolName = RDKitUtil.GetMolName(Mol, MolCount)
     117             MiscUtil.PrintWarning("Ignoring empty molecule: %s" % MolName)
     118             continue
     119
     120         ValidMolCount += 1
     121         MolWithHs = Chem.AddHs(Mol)
     122
     123         # Retrieve charges...
     124         PartialCharges = CalculateMolPartialCharges(MolWithHs, MolCount)
     125         if not len(PartialCharges):
     126             CalcFailedCount += 1
     127             continue
     128
     129         # Write out charges...
     130         WriteMolPartialCharges(Writer, MolWithHs, PartialCharges, Compute2DCoords)
     131
     132     if Writer is not None:
     133         Writer.close()
     134
     135     MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount)
     136     MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount)
     137     MiscUtil.PrintInfo("Number of molecules failed during calculation of partial charges: %d" % CalcFailedCount)
     138     MiscUtil.PrintInfo("Number of ignored molecules: %d" % (MolCount - ValidMolCount + CalcFailedCount))
     139
     140 def CalculateMolPartialCharges(Mol, MolCount):
     141     """Calculate partial atomic charges for a molecule."""
     142
     143     PartialCharges = []
     144     if OptionsInfo["MMFFChargesMode"]:
     145         if  AllChem.MMFFHasAllMoleculeParams(Mol):
     146             MMFFProp = AllChem.MMFFGetMoleculeProperties(Mol)
     147             PartialCharges = [MMFFProp.GetMMFFPartialCharge(AtomIndex) for AtomIndex in range(Mol.GetNumAtoms())]
     148         else:
     149             MolName = RDKitUtil.GetMolName(Mol, MolCount)
     150             MiscUtil.PrintWarning("Failed to calculate MMFF partial charges for molecule, %s: Missing forcefield parameters" % MolName)
     151             return PartialCharges
     152     else:
     153         rdPartialCharges.ComputeGasteigerCharges(Mol, nIter = OptionsInfo["NumIters"], throwOnParamFailure = OptionsInfo["AllowParamFailure"])
     154         PartialCharges = [Atom.GetProp("_GasteigerCharge") for Atom in Mol.GetAtoms()]
     155
     156     # Format charges...
     157     PartialCharges = ["%.*f" % (OptionsInfo["Precision"], float(Value)) for Value in PartialCharges]
     158
     159     return PartialCharges
     160
     161 def WriteMolPartialCharges(Writer, Mol, PartialCharges, Compute2DCoords):
     162     """Write out partial atomic charges for a molecule."""
     163
     164     if OptionsInfo["AtomAliasesFormatMode"]:
     165         for Atom, PartialCharge in zip(Mol.GetAtoms(), PartialCharges):
     166             Atom.SetProp('molFileAlias', PartialCharge)
     167     else:
     168         ChargesValues = "\n".join(PartialCharges)
     169         Mol.SetProp(OptionsInfo["DataFieldLabel"], ChargesValues)
     170
     171     if Compute2DCoords:
     172         AllChem.Compute2DCoords(Mol)
     173
     174     Writer.write(Mol)
     175
     176 def ProcessOptions():
     177     """Process and validate command line arguments and options"""
     178
     179     MiscUtil.PrintInfo("Processing options...")
     180
     181     # Validate options...
     182     ValidateOptions()
     183
     184     AllowParamFailure = True
     185     if re.match("^No", Options["--allowParamFailure"], re.I):
     186         AllowParamFailure = False
     187     OptionsInfo["AllowParamFailure"] = AllowParamFailure
     188
     189     AtomAliasesFormatMode = True
     190     if re.match("^DataField", Options["--chargesSDFormat"], re.I):
     191         AtomAliasesFormatMode = False
     192     OptionsInfo["AtomAliasesFormatMode"] = AtomAliasesFormatMode
     193
     194     OptionsInfo["DataFieldLabel"] = Options["--dataFieldLabel"]
     195
     196     MMFFChargesMode = False
     197     if re.match("^MMFF", Options["--mode"], re.I):
     198         MMFFChargesMode = True
     199     OptionsInfo["Mode"] = Options["--mode"]
     200     OptionsInfo["MMFFChargesMode"] = MMFFChargesMode
     201
     202     OptionsInfo["Infile"] = Options["--infile"]
     203     OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters("--infileParams", Options["--infileParams"], Options["--infile"])
     204
     205     OptionsInfo["Outfile"] = Options["--outfile"]
     206     OptionsInfo["OutfileParams"] = MiscUtil.ProcessOptionOutfileParameters("--outfileParams", Options["--outfileParams"], Options["--infile"], Options["--outfile"])
     207
     208     OptionsInfo["Overwrite"] = Options["--overwrite"]
     209
     210     OptionsInfo["NumIters"] = int(Options["--numIters"])
     211     OptionsInfo["Precision"] = int(Options["--precision"])
     212
     213 def RetrieveOptions():
     214     """Retrieve command line arguments and options"""
     215
     216     # Get options...
     217     global Options
     218     Options = docopt(_docoptUsage_)
     219
     220     # Set current working directory to the specified directory...
     221     WorkingDir = Options["--workingdir"]
     222     if WorkingDir:
     223         os.chdir(WorkingDir)
     224
     225     # Handle examples option...
     226     if "--examples" in Options and Options["--examples"]:
     227         MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
     228         sys.exit(0)
     229
     230 def ValidateOptions():
     231     """Validate option values"""
     232
     233     MiscUtil.ValidateOptionTextValue("-a, --allowParamFailure", Options["--allowParamFailure"], "yes no")
     234     MiscUtil.ValidateOptionTextValue("-c, --chargesSDFormat", Options["--chargesSDFormat"], "AtomAliases DataField")
     235
     236     MiscUtil.ValidateOptionTextValue("-m, --mode", Options["--mode"], "Gasteiger MMFF")
     237
     238     MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
     239     MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "sdf sd mol smi csv tsv txt")
     240
     241     MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "sdf sd")
     242     MiscUtil.ValidateOptionsOutputFileOverwrite("-o, --outfile", Options["--outfile"], "--overwrite", Options["--overwrite"])
     243     MiscUtil.ValidateOptionsDistinctFileNames("-i, --infile", Options["--infile"], "-o, --outfile", Options["--outfile"])
     244
     245     MiscUtil.ValidateOptionIntegerValue("-n, --numIters", Options["--numIters"], {">": 0})
     246     MiscUtil.ValidateOptionIntegerValue("-p, --precision", Options["--precision"], {">": 0})
     247
     248 # Setup a usage string for docopt...
     249 _docoptUsage_ = """
     250 RDKitCalculatePartialCharges.py - Calculate partial atomic charges
     251
     252 Usage:
     253     RDKitCalculatePartialCharges.py [--allowParamFailure <yes or no>]
     254                                           [--chargesSDFormat <AtomAliases or DataField>]  [--dataFieldLabel <text>]
     255                                           [--infileParams <Name,Value,...>] [--mode <Gasteiger or MMFF>]  [--numIters <number>]
     256                                           [--outfileParams <Name,Value,...>] [--precision <number>] [--overwrite]
     257                                           [-w <dir>] -i <infile> -o <outfile>
     258     RDKitCalculatePartialCharges.py -h | --help | -e | --examples
     259
     260 Description:
     261     Calculate partial charges for atoms in molecules and write them out to a SD file.
     262     The hydrogens are automatically added to molecules before calculating partial
     263     charges.
     264
     265     The supported input file formats are: Mol (.mol), SD (.sdf, .sd), SMILES (.smi,
     266     .txt, .csv, .tsv)
     267
     268     The supported output file format are: SD File (.sdf, .sd)
     269
     270 Options:
     271     -a, --allowParamFailure <yes or no>  [default: yes]
     272         Allow calculation of Gasteiger partial charges to proceed for molecules
     273         containing atoms with unknown parameters. The atoms with unknown
     274         parameters are removed from the calculations by setting their values to
     275         zero.
     276     -c, --chargesSDFormat <AtomAliases or DataField>  [default: AtomAliases]
     277         Format for writing out partial atomic charges to SD file. Possible values:
     278         AtomAliases or DataField.
     279
     280         The charges are stored as atom property named 'molFileAlias' for
     281         'AtomAliases' format and may be retrieved using the RDKit function
     282         'GetProp' for atoms: Aotm.GetProp('molFileAliases').
     283
     284         The charges are stored under a data field label speciefied using
     285         '-d, --dataFieldLabel' for 'DataField' format and may be retrieved using the
     286         RDKit function 'GetProp' for molecues.
     287     -d, --dataFieldLabel <text>  [default: PartialCharges]
     288         Data field label to use for storing charged in SD file during 'DataField' value
     289         of '-c, --chargesSDFormat'.
     290     -e, --examples
     291         Print examples.
     292     -h, --help
     293         Print this help message.
     294     -i, --infile <infile>
     295         Input file name.
     296     --infileParams <Name,Value,...>  [default: auto]
     297         A comma delimited list of parameter name and value pairs for reading
     298         molecules from files. The supported parameter names for different file
     299         formats, along with their default values, are shown below:
     300
     301             SD, MOL: removeHydrogens,yes,sanitize,yes,strictParsing,yes
     302             SMILES: smilesColumn,1,smilesNameColumn,2,smilesDelimiter,space,
     303                 smilesTitleLine,auto,sanitize,yes
     304
     305         Possible values for smilesDelimiter: space, comma or tab.
     306     -m, --mode <Gasteiger or MMFF>  [default: Gasteiger]
     307         Type of partial atomic charges to calculate. Possible values: Gasteiger
     308         [ Ref 138 ] or Merk Molecular Mechanics Fore Field (MMFF) [ Ref 83-87 ].
     309     -n, --numIters <number>  [default: 12]
     310         Number of iterations to perform during calculation of Gasteiger charges.
     311     -o, --outfile <outfile>
     312         Output file name.
     313     --outfileParams <Name,Value,...>  [default: auto]
     314         A comma delimited list of parameter name and value pairs for writing
     315         molecules to files. The supported parameter names for different file
     316         formats, along with their default values, are shown below:
     317
     318             SD: compute2DCoords,auto,kekulize,no
     319
     320         Default value for compute2DCoords: yes for SMILES input file; no for all other
     321         file types.
     322     -p, --precision <number>  [default: 3]
     323         Floating point precision for writing the calculated partial atomic charges.
     324     --overwrite
     325         Overwrite existing files.
     326     -w, --workingdir <dir>
     327         Location of working directory which defaults to the current directory.
     328
     329 Examples:
     330     To calculate Gasteiger partial atomic charges for molecules in a SMILES
     331     file and write them out to a SD file as atom aliases, type:
     332
     333         % RDKitCalculatePartialCharges.py  -i Sample.smi -o SampleOut.sdf
     334
     335     To calculate MMFF forcefield partial atomic charges for molecules in a SD
     336     file and write them out to a SD file under 'PartialCharges' data field, type:
     337
     338         % RDKitCalculatePartialCharges.py  -m MMFF -c DataField -i Sample.sdf
     339           -o SampleOut.sdf
     340
     341     To calculate Gasteiger partial atomic charges for molecules in a SMILES
     342     file and write them out to a SD file under a data field named 'GasteigerCharges',
     343     type:
     344
     345         % RDKitCalculatePartialCharges.py  -m Gasteiger -c DataField
     346           -d GasteigerCharges -p 4 -i Sample.smi -o SampleOut.sdf
     347
     348     To calculate Gasteiger partial atomic charges for molecules in a CSV SMILES
     349     file, SMILES strings in column 1, name in column 2, and write out a SD file
     350     containing charges as atom aliases, type:
     351
     352         % RDKitCalculatePartialCharges.py --infileParams
     353           "smilesDelimiter,comma,smilesTitleLine,yes,smilesColumn,1,
     354           smilesNameColumn,2" --outfileParams "compute2DCoords,yes"
     355           -i SampleSMILES.csv -o SampleOut.sdf
     356
     357 Author:
     358     Manish Sud(msud@san.rr.com)
     359
     360 See also:
     361     RDKitCalculateMolecularDescriptors.py, RDKitCalculateRMSD.py,
     362     RDKitCompareMoleculeShapes.py, RDKitConvertFileFormat.py,
     363
     364 Copyright:
     365     Copyright (C) 2018 Manish Sud. All rights reserved.
     366
     367     The functionality available in this script is implemented using RDKit, an
     368     open source toolkit for cheminformatics developed by Greg Landrum.
     369
     370     This file is part of MayaChemTools.
     371
     372     MayaChemTools is free software; you can redistribute it and/or modify it under
     373     the terms of the GNU Lesser General Public License as published by the Free
     374     Software Foundation; either version 3 of the License, or (at your option) any
     375     later version.
     376
     377 """
     378
     379 if __name__ == "__main__":
     380     main()