       1 #!/bin/env python
       2 #
       3 # File: RDKitSearchSMARTS.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 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     PerformSearch()
      76
      77     MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
      78     MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
      79
      80 def PerformSearch():
      81     """Perform search using specified SMARTS pattern."""
      82
      83     Infile = OptionsInfo["Infile"]
      84     Outfile = OptionsInfo["Outfile"]
      85
      86     CountMode = OptionsInfo["CountMode"]
      87     NegateMatch = OptionsInfo["NegateMatch"]
      88     UseChirality = OptionsInfo["UseChirality"]
      89
      90     # Set up a pattern molecule...
      91     PatternMol = Chem.MolFromSmarts(OptionsInfo["Pattern"])
      92
      93     # Setup a molecule reader...
      94     MiscUtil.PrintInfo("\nProcessing file %s..." % Infile)
      95     Mols  = RDKitUtil.ReadMolecules(Infile, **OptionsInfo["InfileParams"])
      96
      97     # Set up a molecule writer...
      98     Writer = None
      99     if not CountMode:
     100         Writer = RDKitUtil.MoleculesWriter(Outfile, **OptionsInfo["OutfileParams"])
     101         if Writer is None:
     102             MiscUtil.PrintError("Failed to setup a writer for output fie %s " % Outfile)
     103         MiscUtil.PrintInfo("Generating file %s..." % Outfile)
     104
     105     # Process molecules...
     106     MolCount = 0
     107     ValidMolCount = 0
     108     MatchCount = 0
     109
     110     Compute2DCoords = OptionsInfo["OutfileParams"]["Compute2DCoords"]
     111
     112     for Mol in Mols:
     113         MolCount += 1
     114
     115         if Mol is None:
     116             continue
     117
     118         if RDKitUtil.IsMolEmpty(Mol):
     119             MolName = RDKitUtil.GetMolName(Mol, MolCount)
     120             MiscUtil.PrintWarning("Ignoring empty molecule: %s" % MolName)
     121             continue
     122
     123         ValidMolCount += 1
     124
     125         MolMatched = Mol.HasSubstructMatch(PatternMol, useChirality = UseChirality)
     126         if MolMatched != NegateMatch:
     127             MatchCount += 1
     128             if not CountMode:
     129                 if Compute2DCoords:
     130                     AllChem.Compute2DCoords(Mol)
     131                 Writer.write(Mol)
     132
     133     if Writer is not None:
     134         Writer.close()
     135
     136     MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount)
     137     MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount)
     138     MiscUtil.PrintInfo("Number of ignored molecules: %d" % (MolCount - ValidMolCount))
     139
     140     if NegateMatch:
     141         MiscUtil.PrintInfo("\nTotal number of not matched molecules: %d" % MatchCount)
     142     else:
     143         MiscUtil.PrintInfo("\nTotal number of matched molecules: %d" % MatchCount)
     144
     145 def ProcessOptions():
     146     """Process and validate command line arguments and options"""
     147
     148     MiscUtil.PrintInfo("Processing options...")
     149
     150     # Validate options...
     151     ValidateOptions()
     152
     153     OptionsInfo["Infile"] = Options["--infile"]
     154     OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters("--infileParams", Options["--infileParams"], Options["--infile"])
     155
     156     OptionsInfo["Outfile"] = Options["--outfile"]
     157     OptionsInfo["OutfileParams"] = MiscUtil.ProcessOptionOutfileParameters("--outfileParams", Options["--outfileParams"], Options["--infile"], Options["--outfile"])
     158
     159     OptionsInfo["Overwrite"] = Options["--overwrite"]
     160
     161     OptionsInfo["CountMode"] = False
     162     if re.match("^count$", Options["--mode"], re.I):
     163         OptionsInfo["CountMode"] = True
     164
     165     OptionsInfo["NegateMatch"] = False
     166     if re.match("^yes$", Options["--negate"], re.I):
     167         OptionsInfo["NegateMatch"] = True
     168
     169     OptionsInfo["Pattern"] = Options["--pattern"]
     170
     171     OptionsInfo["UseChirality"] = False
     172     if re.match("^yes$", Options["--useChirality"], re.I):
     173         OptionsInfo["UseChirality"] = True
     174
     175 def RetrieveOptions():
     176     """Retrieve command line arguments and options"""
     177
     178     # Get options...
     179     global Options
     180     Options = docopt(_docoptUsage_)
     181
     182     # Set current working directory to the specified directory...
     183     WorkingDir = Options["--workingdir"]
     184     if WorkingDir:
     185         os.chdir(WorkingDir)
     186
     187     # Handle examples option...
     188     if "--examples" in Options and Options["--examples"]:
     189         MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
     190         sys.exit(0)
     191
     192 def ValidateOptions():
     193     """Validate option values"""
     194
     195     MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
     196     MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "sdf sd smi txt csv tsv")
     197     if Options["--outfile"]:
     198         MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "sdf sd smi")
     199         MiscUtil.ValidateOptionsOutputFileOverwrite("-o, --outfile", Options["--outfile"], "--overwrite", Options["--overwrite"])
     200         MiscUtil.ValidateOptionsDistinctFileNames("-i, --infile", Options["--infile"], "-o, --outfile", Options["--outfile"])
     201
     202     MiscUtil.ValidateOptionTextValue("-m, --mode", Options["--mode"], "retrieve count")
     203     if re.match("^retrieve$", Options["--mode"], re.I):
     204         if not Options["--outfile"]:
     205             MiscUtil.PrintError("The outfile must be specified using \"-o, --outfile\" during \"retrieve\" value of \"-m, --mode\" option")
     206
     207     MiscUtil.ValidateOptionTextValue("-n, --negate", Options["--negate"], "yes no")
     208
     209     PatternMol = Chem.MolFromSmarts(Options["--pattern"])
     210     if PatternMol is None:
     211         MiscUtil.PrintError("The value specified, %s, using option \"-p, --pattern\" is not a valid SMARTS: Failed to create pattern molecule" % Options["--pattern"])
     212
     213     MiscUtil.ValidateOptionTextValue("--useChirality", Options["--useChirality"], "yes no")
     214
     215 # Setup a usage string for docopt...
     216 _docoptUsage_ = """
     217 RDKitSearchSMARTS.py - Perform a substructure search using SMARTS pattern
     218
     219 Usage:
     220     RDKitSearchSMARTS.py  [--infileParams <Name,Value,...>] [--mode <retrieve or count>]
     221                           [--negate <yes or no>] [--outfileParams <Name,Value,...>] [--overwrite]
     222                           [--useChirality <yes or no>] [-w <dir>] [-o <outfile>] -p <SMARTS> -i <infile>
     223     RDKitSearchSMARTS.py -h | --help | -e | --examples
     224
     225 Description:
     226     Perform a substructure search in an input file using specified SMARTS pattern and
     227     write out the matched molecules to an output file or simply count the number
     228     of matches.
     229
     230     The supported input file formats are: SD (.sdf, .sd), SMILES (.smi., csv, .tsv, .txt)
     231
     232     The supported output file formats are: SD (.sdf, .sd), SMILES (.smi)
     233
     234 Options:
     235     -e, --examples
     236         Print examples.
     237     -h, --help
     238         Print this help message.
     239     -i, --infile <infile>
     240         Input file name.
     241     --infileParams <Name,Value,...>  [default: auto]
     242         A comma delimited list of parameter name and value pairs for reading
     243         molecules from files. The supported parameter names for different file
     244         formats, along with their default values, are shown below:
     245
     246             SD, MOL: removeHydrogens,yes,sanitize,yes,strictParsing,yes
     247             SMILES: smilesColumn,1,smilesNameColumn,2,smilesDelimiter,space,
     248                 smilesTitleLine,auto,sanitize,yes
     249
     250         Possible values for smilesDelimiter: space, comma or tab.
     251     -m, --mode <retrieve or count>  [default: retrieve]
     252         Specify whether to retrieve and write out matched molecules to an output
     253         file or simply count the number of matches.
     254     -n, --negate <yes or no>  [default: no]
     255         Specify whether to find molecules not matching the specified SMARTS pattern.
     256     -o, --outfile <outfile>
     257         Output file name.
     258     --outfileParams <Name,Value,...>  [default: auto]
     259         A comma delimited list of parameter name and value pairs for writing
     260         molecules to files. The supported parameter names for different file
     261         formats, along with their default values, are shown below:
     262
     263             SD: compute2DCoords,auto,kekulize,no
     264             SMILES: kekulize,no,smilesDelimiter,space, smilesIsomeric,yes,
     265                 smilesTitleLine,yes
     266
     267         Default value for compute2DCoords: yes for SMILES input file; no for all other
     268         file types.
     269     --overwrite
     270         Overwrite existing files.
     271     -p, --pattern <SMARTS>  [default: none]
     272         SMARTS pattern for performing search.
     273     -u, --useChirality <yes or no>  [default: no]
     274         Use stereochemistry information for SMARTS search.
     275     -w, --workingdir <dir>
     276         Location of working directory which defaults to the current directory.
     277
     278 Examples:
     279     To retrieve molecules containing the substructure corresponding to a specified
     280     SMARTS pattern and write out a SMILES file, type:
     281
     282         % RDKitSearchSMARTS.py -p 'c1ccccc1' -i Sample.smi -o SampleOut.smi
     283
     284     To only count the number of molecules containing the substructure corresponding
     285     to a specified SMARTS pattern without writing out any file, type:
     286
     287         % RDKitSearchSMARTS.py -m count -p 'c1ccccc1' -i Sample.smi
     288
     289     To count the number of molecules in a SD file not containing the substructure
     290     corresponding to a specified SMARTS pattern and write out a SD file, type:
     291
     292         % RDKitSearchSMARTS.py -n yes -p 'c1ccccc1' -i Sample.sdf -o SampleOut.sdf
     293
     294     To retrieve molecules containing the substructure corresponding to a specified
     295     SMARTS pattern from a CSV SMILES file, SMILES strings in column 1, name in
     296     and write out a SD file, type:
     297
     298         % RDKitSearchSMARTS.py -p 'c1ccccc1' --infileParams
     299           "smilesDelimiter,comma,smilesTitleLine,yes,smilesColumn,1,
     300           smilesNameColumn,2" --outfileParams "compute2DCoords,yes"
     301           -i SampleSMILES.csv -o SampleOut.sdf
     302
     303 Author:
     304     Manish Sud(msud@san.rr.com)
     305
     306 See also:
     307     RDKitConvertFileFormat.py, RDKitFilterPAINS.py, RDKitSearchFunctionalGroups.py
     308
     309 Copyright:
     310     Copyright (C) 2018 Manish Sud. All rights reserved.
     311
     312     The functionality available in this script is implemented using RDKit, an
     313     open source toolkit for cheminformatics developed by Greg Landrum.
     314
     315     This file is part of MayaChemTools.
     316
     317     MayaChemTools is free software; you can redistribute it and/or modify it under
     318     the terms of the GNU Lesser General Public License as published by the Free
     319     Software Foundation; either version 3 of the License, or (at your option) any
     320     later version.
     321
     322 """
     323
     324 if __name__ == "__main__":
     325     main()