MayaChemTools

   1 #!/bin/env python
   2 #
   3 # File: RDKitCalculateRMSD.py
   4 # Author: Manish Sud <msud@san.rr.com>
   5 #
   6 # Copyright (C) 2019 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, rdMolAlign
  42 except ImportError as ErrMsg:
  43     sys.stderr.write("\nFailed to import RDKit module/package: %s\n" % ErrMsg)
  44     sys.stderr.write("Check/update your RDKit environment and try again.\n\n")
  45     sys.exit(1)
  46 
  47 # MayaChemTools imports...
  48 try:
  49     from docopt import docopt
  50     import MiscUtil
  51     import RDKitUtil
  52 except ImportError as ErrMsg:
  53     sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg)
  54     sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n")
  55     sys.exit(1)
  56 
  57 ScriptName = os.path.basename(sys.argv[0])
  58 Options = {}
  59 OptionsInfo = {}
  60 
  61 def main():
  62     """Start execution of the script"""
  63     
  64     MiscUtil.PrintInfo("\n%s (RDK v%s; %s): Starting...\n" % (ScriptName, rdBase.rdkitVersion, time.asctime()))
  65     
  66     (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime()
  67     
  68     # Retrieve command line arguments and options...
  69     RetrieveOptions()
  70     
  71     # Process and validate command line arguments and options...
  72     ProcessOptions()
  73     
  74     # Perform actions required by the script...
  75     CalculateRMSD()
  76     
  77     MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
  78     MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
  79 
  80 def CalculateRMSD():
  81     """Calculate RMSD values."""
  82     
  83     if not re.match("^(OneToOne|AllToAll|FirstToAll)$", OptionsInfo["Mode"], re.I):
  84         MiscUtil.PrintError("RMSD couldn't be calculated: Specified mode, %s, is not supported" % OptionsInfo["Mode"])
  85         
  86     RefFile = OptionsInfo["RefFile"]
  87     ProbeFile = OptionsInfo["ProbeFile"]
  88     
  89     Outfile = OptionsInfo["Outfile"]
  90     OutDelim = OptionsInfo["OutDelim"]
  91 
  92     # Read reference and probe molecules...
  93     OptionsInfo["InfileParams"]["AllowEmptyMols"] = False
  94     
  95     MiscUtil.PrintInfo("\nProcessing file %s..." % (RefFile))
  96     ValidRefMols, RefMolCount, ValidRefMolCount  = RDKitUtil.ReadAndValidateMolecules(RefFile, **OptionsInfo["InfileParams"])
  97     
  98     MiscUtil.PrintInfo("Processing file %s..." % (ProbeFile))
  99     ValidProbeMols, ProbeMolCount, ValidProbeMolCount  = RDKitUtil.ReadAndValidateMolecules(ProbeFile, **OptionsInfo["InfileParams"])
 100 
 101     # Set up output file...
 102     MiscUtil.PrintInfo("\nGenerating file %s...\n" % Outfile)
 103     OutFH = open(Outfile, "w")
 104     if OutFH is None:
 105         MiscUtil.PrintError("Couldn't open output file: %s.\n" % (Outfile))
 106     
 107     Line = "RefMolID%sProbeMolID%sRMSD\n" % (OutDelim, OutDelim)
 108     OutFH.write(Line)
 109     
 110     if re.match("^OneToOne$", OptionsInfo["Mode"], re.I):
 111         CalculateOneToOneRMSDValues(ValidRefMols, ValidProbeMols, OutFH, OutDelim)
 112     elif re.match("^AllToAll$", OptionsInfo["Mode"], re.I):
 113         CalculateAllToAllRMSDValues(ValidRefMols, ValidProbeMols, OutFH, OutDelim)
 114     elif re.match("^FirstToAll$", OptionsInfo["Mode"], re.I):
 115         CalculateFirstToAllRMSDValues(ValidRefMols, ValidProbeMols, OutFH, OutDelim)
 116     else:
 117         MiscUtil.PrintError("RMSD couldn't be calculated: Specified mode, %s, is not supported" % OptionsInfo["Mode"])
 118 
 119     OutFH.close()
 120     
 121     MiscUtil.PrintInfo("\nTotal number of molecules: Reference - %d; Probe - %d" % (RefMolCount, ProbeMolCount))
 122     MiscUtil.PrintInfo("Number of valid molecules: Reference - %d; Probe - %d" % (ValidRefMolCount, ValidProbeMolCount))
 123     MiscUtil.PrintInfo("Number of ignored molecules:  Reference - %d; Probe - %d" % ((RefMolCount - ValidRefMolCount), (ProbeMolCount - ValidProbeMolCount)))
 124     
 125 def CalculateOneToOneRMSDValues(ValidRefMols, ValidProbeMols, OutFH, OutDelim):
 126     """Calculate pairwise RMSD values."""
 127     
 128     ValidRefMolCount = len(ValidRefMols)
 129     ValidProbeMolCount = len(ValidProbeMols)
 130     
 131     MolCount = ValidRefMolCount
 132     if ValidRefMolCount > ValidProbeMolCount:
 133         MolCount = ValidProbeMolCount
 134     
 135     if ValidRefMolCount != ValidProbeMolCount:
 136         MiscUtil.PrintWarning("Number of valid reference molecules, %d,  is not equal to number of valid probe molecules, %d .\n" % (ValidRefMolCount, ValidProbeMolCount))
 137         MiscUtil.PrintWarning("Pairwise RMSD will be calculated only for first %s molecules.\n" % (MolCount))
 138 
 139     # Process molecules...
 140     for MolIndex in range(0, MolCount):
 141         RefMol = ValidRefMols[MolIndex]
 142         ProbeMol = ValidProbeMols[MolIndex]
 143 
 144         RefMolName = RDKitUtil.GetMolName(RefMol, (MolIndex + 1))
 145         ProbeMolName = RDKitUtil.GetMolName(ProbeMol, (MolIndex + 1))
 146 
 147         RMSD = CalculateRMSDValue(RefMol, ProbeMol)
 148         
 149         Line = "%s%s%s%s%s\n" % (RefMolName, OutDelim, ProbeMolName, OutDelim, RMSD)
 150         OutFH.write(Line)
 151         
 152 def CalculateAllToAllRMSDValues(ValidRefMols, ValidProbeMols, OutFH, OutDelim):
 153     """Calculate RMSD values between all pairs of molecules."""
 154     
 155     # Process molecules...
 156     RefMolCount = 0
 157     for RefMol in ValidRefMols:
 158         RefMolCount += 1
 159         RefMolName = RDKitUtil.GetMolName(RefMol, RefMolCount)
 160 
 161         ProbeMolCount = 0
 162         for ProbeMol in ValidProbeMols:
 163             ProbeMolCount += 1
 164             ProbeMolName = RDKitUtil.GetMolName(ProbeMol, ProbeMolCount)
 165             
 166             RMSD = CalculateRMSDValue(RefMol, ProbeMol)
 167 
 168             Line = "%s%s%s%s%s\n" % (RefMolName, OutDelim, ProbeMolName, OutDelim, RMSD)
 169             OutFH.write(Line)
 170         
 171 def CalculateFirstToAllRMSDValues(ValidRefMols, ValidProbeMols, OutFH, OutDelim):
 172     """Calculate RMSD values between first reference molecues and all probe molecules. """
 173     
 174     # Process molecules...
 175     RefMol = ValidRefMols[0]
 176     RefMolCount = 1
 177     RefMolName = RDKitUtil.GetMolName(RefMol, RefMolCount)
 178 
 179     ProbeMolCount = 0
 180     for ProbeMol in ValidProbeMols:
 181         ProbeMolCount += 1
 182         ProbeMolName = RDKitUtil.GetMolName(ProbeMol, ProbeMolCount)
 183             
 184         RMSD = CalculateRMSDValue(RefMol, ProbeMol)
 185 
 186         Line = "%s%s%s%s%s\n" % (RefMolName, OutDelim, ProbeMolName, OutDelim, RMSD)
 187         OutFH.write(Line)
 188         
 189 def CalculateRMSDValue(RefMol, ProbeMol):
 190     """Calculate RMSD value for a pair of molecules and return it as a string"""
 191 
 192     try:
 193         if OptionsInfo["UseBestRMSD"]:
 194             RMSD = AllChem.GetBestRMS(RefMol, ProbeMol)
 195         else:
 196             RMSD = rdMolAlign.AlignMol(ProbeMol, RefMol, maxIters = OptionsInfo["MaxIters"])
 197         RMSD = "%.2f" % RMSD
 198     except (RuntimeError, ValueError):
 199         RMSD = "None"
 200         
 201     return RMSD
 202 
 203 def ProcessOptions():
 204     """Process and validate command line arguments and options"""
 205     
 206     MiscUtil.PrintInfo("Processing options...")
 207     
 208     # Validate options...
 209     ValidateOptions()
 210     
 211     OptionsInfo["CalcRMSD"] = Options["--calcRMSD"]
 212     OptionsInfo["UseBestRMSD"] = False
 213     if re.match("^BestRMSD$", OptionsInfo["CalcRMSD"], re.I):
 214         OptionsInfo["UseBestRMSD"] = True
 215     
 216     OptionsInfo["MaxIters"] = int(Options["--maxIters"])
 217     
 218     OptionsInfo["Mode"] = Options["--mode"]
 219     
 220     OptionsInfo["RefFile"] = Options["--reffile"]
 221     OptionsInfo["ProbeFile"] = Options["--probefile"]
 222 
 223     # No need for any RDKit specific --outfileParams....
 224     OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters("--infileParams", Options["--infileParams"])
 225     
 226     OptionsInfo["Outfile"] = Options["--outfile"]
 227     
 228     OptionsInfo["Overwrite"] = Options["--overwrite"]
 229     
 230     OptionsInfo["OutDelim"] = " "
 231     if MiscUtil.CheckFileExt(OptionsInfo["Outfile"], "csv"):
 232         OptionsInfo["OutDelim"] = ","
 233     elif MiscUtil.CheckFileExt(OptionsInfo["Outfile"], "tsv txt"):
 234         OptionsInfo["OutDelim"] = "\t"
 235     else:
 236         MiscUtil.PrintError("The file name specified , %s, for option \"--outfile\" is not valid. Supported file formats: csv tsv txt\n" % (OptionsInfo["Outfile"]))
 237 
 238 def RetrieveOptions():
 239     """Retrieve command line arguments and options"""
 240     
 241     # Get options...
 242     global Options
 243     Options = docopt(_docoptUsage_)
 244     
 245     # Set current working directory to the specified directory...
 246     WorkingDir = Options["--workingdir"]
 247     if WorkingDir:
 248         os.chdir(WorkingDir)
 249     
 250     # Handle examples option...
 251     if "--examples" in Options and Options["--examples"]:
 252         MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
 253         sys.exit(0)
 254 
 255 def ValidateOptions():
 256     """Validate option values"""
 257     
 258     MiscUtil.ValidateOptionTextValue("--calcRMSD", Options["--calcRMSD"], "RMSD BestRMSD")
 259     
 260     MiscUtil.ValidateOptionIntegerValue("--maxIters", Options["--maxIters"], {">": 0})
 261     MiscUtil.ValidateOptionTextValue("--mode", Options["--mode"], "OneToOne  AllToAll FirstToAll")
 262     
 263     MiscUtil.ValidateOptionFilePath("-r, --reffile", Options["--reffile"])
 264     MiscUtil.ValidateOptionFileExt("-r, --reffile", Options["--reffile"], "sdf sd mol")
 265     
 266     MiscUtil.ValidateOptionFilePath("-p, --probefile", Options["--probefile"])
 267     MiscUtil.ValidateOptionFileExt("-p, --probefile", Options["--probefile"], "sdf sd mol")
 268         
 269     MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "csv tsv txt")
 270     MiscUtil.ValidateOptionsOutputFileOverwrite("-o, --outfile", Options["--outfile"], "--overwrite", Options["--overwrite"])
 271     MiscUtil.ValidateOptionsDistinctFileNames("-r, --reffile", Options["--reffile"], "-o, --outfile", Options["--outfile"])
 272     MiscUtil.ValidateOptionsDistinctFileNames("-p, --probefile", Options["--probefile"], "-o, --outfile", Options["--outfile"])
 273     
 274 # Setup a usage string for docopt...
 275 _docoptUsage_ = """
 276 RDKitCalculateRMSD.py - Calculate RMSD between molecules
 277 
 278 Usage:
 279     RDKitCalculateRMSD.py [--calcRMSD <RMSD, BestRMSD>] [--infileParams <Name,Value,...>]
 280                           [--maxIters <number>] [--mode <OneToOne, AllToAll, FirstToAll>]
 281                           [--overwrite] [-w <dir>] -r <reffile> -p <probefile> -o <outfile> 
 282     RDKitCalculateRMSD.py -h | --help | -e | --examples
 283 
 284 Description:
 285     Calculate Root Mean Square Distance (RMSD) between a set of similar molecules in
 286     reference and probe input files. The RDKit function fails to calculate RMSD values for
 287     dissimilar molecules. Consequently, a text string 'None' is written out as a RMSD value
 288     for dissimilar molecule pairs.
 289 
 290     The supported input file formats are: Mol (.mol), SD (.sdf, .sd)
 291 
 292     The supported output file formats are:  CSV (.csv), TSV (.tsv, .txt)
 293 
 294 Options:
 295     -c, --calcRMSD <RMSD, BestRMSD>  [default: RMSD]
 296         Methodology for calculating RMSD values. Possible values: RMSD, BestRMSD.
 297         During BestRMSMode mode, the RDKit 'function AllChem.GetBestRMS' is used to
 298         align and calculate RMSD. This function calculates optimal RMSD for aligning two
 299         molecules, taking symmetry into account. Otherwise, the RMSD value is calculated
 300         using 'AllChem.AlignMol function' without changing the atom order. A word to the
 301         wise from RDKit documentation: The AllChem.GetBestRMS function will attempt to
 302         align all permutations of matching atom orders in both molecules, for some molecules
 303         it will lead to 'combinatorial explosion'.
 304     --infileParams <Name,Value,...>  [default: auto]
 305         A comma delimited list of parameter name and value pairs for reading
 306         molecules from files. The supported parameter names for different file
 307         formats, along with their default values, are shown below:
 308             
 309             SD, MOL: removeHydrogens,yes,sanitize,yes,strictParsing,yes
 310             
 311     --maxIters <number>  [default: 50]
 312         Maximum number of iterations to perform for each molecule pair during minimization
 313         of RMSD values. This option is ignored during BestRMSD mode.
 314     -m, --mode <OneToOne, AllToAll, FirstToAll>  [default: OneToOne]
 315         Specify how molecules are handled in reference and probe input files during
 316         calculation of RMSD between reference and probe molecules.  Possible values:
 317         OneToOne, AllToAll and AllToFirst. For OneToOne mode, the number of molecules
 318         in reference file must be equal to the number of molecules in probe file. The RMSD
 319         is calculated for each pair of molecules in the reference and probe file and written
 320         to the output file. For AllToAll mode, the RMSD is calculated for each reference
 321         molecule against all probe molecules. For FirstToAll mode, however, the RMSD
 322         is only calculated for the first reference molecule against all probe molecules.
 323     -e, --examples
 324         Print examples.
 325     -h, --help
 326         Print this help message.
 327     -p, --probefile <probefile>
 328         Probe input file name.
 329     -r, --reffile <reffile>
 330         Reference input file name.
 331     -o, --outfile <outfile>
 332         Output file name for writing out RMSD values. Supported text file extensions: csv or tsv.
 333     --overwrite
 334         Overwrite existing files.
 335     -w, --workingdir <dir>
 336         Location of working directory which defaults to the current directory.
 337 
 338 Examples:
 339     To calculate RMSD between pair of molecules in reference and probe molecules in
 340     3D SD files and write out a CSV file containing calculated RMSD values along with
 341     appropriate molecule IDs, type:
 342 
 343         % RDKitCalculateRMSD.py  -r Sample3DRef.sdf -p Sample3DProb.sdf
 344           -o SampleOut.csv
 345 
 346     To calculate RMSD between all molecules in reference and probe molecules in
 347     3D SD files and write out a CSV file containing calculated RMSD values along with
 348     appropriate molecule IDs, type:
 349 
 350         % RDKitCalculateRMSD.py  -m AllToAll -r Sample3DRef.sdf -p
 351           Sample3DProb.sdf -o SampleOut.csv
 352 
 353     To calculate best RMSD between first  molecule in reference all probe molecules
 354     in 3D SD files and write out a TSV file containing calculated RMSD values along with
 355     appropriate molecule IDs, type:
 356 
 357         % RDKitCalculateRMSD.py  -m FirstToAll --calcRMSD BestRMSD -r
 358           Sample3DRef.sdf -p Sample3DProb.sdf -o SampleOut.tsv
 359 
 360     To calculate RMSD between all molecules in reference and probe molecules in
 361     3D SD files without removing hydrogens and write out a TSV file containing
 362     calculated RMSD values along with appropriate molecule IDs, type:
 363 
 364         % RDKitCalculateRMSD.py  -m AllToAll --infileParams
 365           "removeHydrogens,no" -r Sample3DRef.sdf  -p Sample3DProb.sdf
 366           -o SampleOut.tsv
 367 
 368 Author:
 369     Manish Sud(msud@san.rr.com)
 370 
 371 See also:
 372     RDKitCalculateMolecularDescriptors.py, RDKitCompareMoleculeShapes.py, RDKitConvertFileFormat.py,
 373     RDKitGenerateConformers.py, RDKitPerformMinimization.py
 374 
 375 Copyright:
 376     Copyright (C) 2019 Manish Sud. All rights reserved.
 377 
 378     The functionality available in this script is implemented using RDKit, an
 379     open source toolkit for cheminformatics developed by Greg Landrum.
 380 
 381     This file is part of MayaChemTools.
 382 
 383     MayaChemTools is free software; you can redistribute it and/or modify it under
 384     the terms of the GNU Lesser General Public License as published by the Free
 385     Software Foundation; either version 3 of the License, or (at your option) any
 386     later version.
 387 
 388 """
 389 
 390 if __name__ == "__main__":
 391     main()