       1 #!/bin/env python
       2 #
       3 # File: RDKitFilterPAINS.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 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 # MayaChemTools imports...
      49 try:
      50     from docopt import docopt
      51     import MiscUtil
      52     import RDKitUtil
      53 except ImportError as ErrMsg:
      54     sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg)
      55     sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n")
      56     sys.exit(1)
      57
      58 ScriptName = os.path.basename(sys.argv[0])
      59 Options = {}
      60 OptionsInfo = {}
      61
      62 def main():
      63     """Start execution of the script"""
      64
      65     MiscUtil.PrintInfo("\n%s (RDK v%s; %s): Starting...\n" % (ScriptName, rdBase.rdkitVersion, time.asctime()))
      66
      67     (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime()
      68
      69     # Retrieve command line arguments and options...
      70     RetrieveOptions()
      71
      72     # Process and validate command line arguments and options...
      73     ProcessOptions()
      74
      75     # Perform actions required by the script...
      76     PerformFiltering()
      77
      78     MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
      79     MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
      80
      81 def PerformFiltering():
      82     """Filter molecules using SMARTS specified in PAINS filter file."""
      83
      84     Infile = OptionsInfo["Infile"]
      85     Outfile = OptionsInfo["Outfile"]
      86
      87     CountMode = OptionsInfo["CountMode"]
      88     NegateMatch = OptionsInfo["NegateMatch"]
      89     PAINSMode = OptionsInfo["PAINSMode"]
      90
      91     # Setup PAINS patterns and pattern mols...
      92     PAINSPatterns = RetrievePAINSPatterns(PAINSMode)
      93     PAINSPatternMols = SetupPAINSPatternMols(PAINSPatterns)
      94
      95     # Setup a molecule reader...
      96     MiscUtil.PrintInfo("\nProcessing file %s..." % Infile)
      97     Mols  = RDKitUtil.ReadMolecules(Infile, **OptionsInfo["InfileParams"])
      98
      99     # Set up a molecule writer...
     100     Writer = None
     101     if not CountMode:
     102         Writer = RDKitUtil.MoleculesWriter(Outfile, **OptionsInfo["OutfileParams"])
     103         if Writer is None:
     104             MiscUtil.PrintError("Failed to setup a writer for output fie %s " % Outfile)
     105         MiscUtil.PrintInfo("Generating file %s..." % Outfile)
     106
     107     Compute2DCoords = OptionsInfo["OutfileParams"]["Compute2DCoords"]
     108
     109     # Process molecules...
     110     MolCount = 0
     111     ValidMolCount = 0
     112     FilteredCount = 0
     113
     114     for Mol in Mols:
     115         MolCount += 1
     116
     117         if Mol is None:
     118             continue
     119
     120         if RDKitUtil.IsMolEmpty(Mol):
     121             MolName = RDKitUtil.GetMolName(Mol, MolCount)
     122             MiscUtil.PrintWarning("Ignoring empty molecule: %s" % MolName)
     123             continue
     124
     125         ValidMolCount += 1
     126
     127         MolMatched = DoesMoleculeContainsPAINSPattern(Mol, PAINSPatternMols)
     128         if MolMatched == NegateMatch:
     129             FilteredCount += 1
     130             if not CountMode:
     131                 if Compute2DCoords:
     132                     AllChem.Compute2DCoords(Mol)
     133                 Writer.write(Mol)
     134
     135     if Writer is not None:
     136         Writer.close()
     137
     138     MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount)
     139     MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount)
     140     MiscUtil.PrintInfo("Number of ignored molecules: %d" % (MolCount - ValidMolCount))
     141
     142     MiscUtil.PrintInfo("\nTotal number of filtered molecules: %d" % FilteredCount)
     143
     144 def DoesMoleculeContainsPAINSPattern(Mol, PAINSPatternMols):
     145     """Check presence of PAINS pattern in the molecule"""
     146
     147     MolMatched = False
     148
     149     for PatternMol in PAINSPatternMols:
     150         if Mol.HasSubstructMatch(PatternMol, useChirality = True):
     151             MolMatched = True
     152             break
     153
     154     return MolMatched
     155
     156 def RetrievePAINSPatterns(PAINSFilterMode):
     157     """Retrieve PAINS patterns for specified PAINS mode"""
     158
     159     MayaChemToolsDataDir = MiscUtil.GetMayaChemToolsLibDataPath()
     160     PAINSFiltersFilePath = os.path.join(MayaChemToolsDataDir, "PAINSFilters.csv")
     161
     162     MiscUtil.PrintInfo("\nRetrieving PAINS SMARTS patterns for PAINS filter type, %s, from file %s" % (PAINSFilterMode, PAINSFiltersFilePath))
     163
     164     if not os.path.exists(PAINSFiltersFilePath):
     165         MiscUtil.PrintError("The PAINS filters file, %s, doesn't exist.\n" % (PAINSFiltersFilePath))
     166
     167     FilterFile = open(PAINSFiltersFilePath, "r")
     168     if FilterFile is None:
     169         MiscUtil.PrintError("Couldn't open PAINS filter file: %s.\n" % (PAINSFiltersFilePath))
     170
     171     # Collect all PAINS filter lines...
     172     HeaderLine = True
     173     FiltersLines = []
     174     for Line in FilterFile:
     175         Line = Line.rstrip()
     176         # Ignore comments...
     177         if re.match("^#", Line, re.I):
     178             continue
     179         # Ignore header line...
     180         if HeaderLine:
     181             HeaderLine = False
     182             continue
     183         FiltersLines.append(Line)
     184
     185     # Process PAINS filter lines using csv reader...
     186     SMARTSPatterns = []
     187
     188     FiltersReader = csv.reader(FiltersLines, delimiter=',', quotechar='"')
     189     for LineWords in FiltersReader:
     190         FilterType = LineWords[0]
     191         ID = LineWords[1]
     192         SMARTS = LineWords[2]
     193
     194         if re.match("^All$", PAINSFilterMode, re.I) or FilterType.lower() == PAINSFilterMode.lower():
     195             SMARTSPatterns.append(SMARTS)
     196
     197     FilterFile.close()
     198
     199     MiscUtil.PrintInfo("Total number of PAINS SMARTS patterns: %d" % (len(SMARTSPatterns)))
     200
     201     return SMARTSPatterns
     202
     203 def SetupPAINSPatternMols(PAINSPatterns):
     204     """Set up PAINS pattern mols for substructure search"""
     205
     206     MiscUtil.PrintInfo("\nSetting up PAINS pattern molecules for performins substructure search...")
     207     PatternMols = []
     208     for Pattern in PAINSPatterns:
     209         PatternMol = Chem.MolFromSmarts(Pattern)
     210         PatternMols.append(PatternMol)
     211
     212     return PatternMols
     213
     214 def ProcessOptions():
     215     """Process and validate command line arguments and options"""
     216
     217     MiscUtil.PrintInfo("Processing options...")
     218
     219     # Validate options...
     220     ValidateOptions()
     221
     222     OptionsInfo["Infile"] = Options["--infile"]
     223     OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters("--infileParams", Options["--infileParams"], Options["--infile"])
     224
     225     OptionsInfo["Outfile"] = Options["--outfile"]
     226     OptionsInfo["OutfileParams"] = MiscUtil.ProcessOptionOutfileParameters("--outfileParams", Options["--outfileParams"], Options["--infile"], Options["--outfile"])
     227
     228     OptionsInfo["Overwrite"] = Options["--overwrite"]
     229
     230     OptionsInfo["CountMode"] = False
     231     if re.match("^count$", Options["--mode"], re.I):
     232         OptionsInfo["CountMode"] = True
     233
     234     OptionsInfo["NegateMatch"] = False
     235     if re.match("^yes$", Options["--negate"], re.I):
     236         OptionsInfo["NegateMatch"] = True
     237
     238     OptionsInfo["PAINSMode"] = Options["--painsMode"]
     239
     240 def RetrieveOptions():
     241     """Retrieve command line arguments and options"""
     242
     243     # Get options...
     244     global Options
     245     Options = docopt(_docoptUsage_)
     246
     247     # Set current working directory to the specified directory...
     248     WorkingDir = Options["--workingdir"]
     249     if WorkingDir:
     250         os.chdir(WorkingDir)
     251
     252     # Handle examples option...
     253     if "--examples" in Options and Options["--examples"]:
     254         MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
     255         sys.exit(0)
     256
     257 def ValidateOptions():
     258     """Validate option values"""
     259
     260     MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
     261     MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "sdf sd smi txt csv tsv")
     262
     263     MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "sdf sd smi")
     264     MiscUtil.ValidateOptionsOutputFileOverwrite("-o, --outfile", Options["--outfile"], "--overwrite", Options["--overwrite"])
     265     MiscUtil.ValidateOptionsDistinctFileNames("-i, --infile", Options["--infile"], "-o, --outfile", Options["--outfile"])
     266
     267     MiscUtil.ValidateOptionTextValue("-m, --mode", Options["--mode"], "filter count")
     268     if re.match("^filter$", Options["--mode"], re.I):
     269         if not Options["--outfile"]:
     270             MiscUtil.PrintError("The outfile must be specified using \"-o, --outfile\" during \"filter\" value of \"-m, --mode\" option")
     271
     272     MiscUtil.ValidateOptionTextValue("-n, --negate", Options["--negate"], "yes no")
     273
     274     MiscUtil.ValidateOptionTextValue("-p, --painsMode", Options["--painsMode"], "All A B C")
     275
     276 # Setup a usage string for docopt...
     277 _docoptUsage_ = """
     278 RDKitFilterPAINS.py - Filter PAINS molecules
     279
     280 Usage:
     281     RDKitFilterPAINS.py  [--infileParams <Name,Value,...>] [--mode <filter or count>]
     282                          [ --outfileParams <Name,Value,...> ] [--painsMode <All, A, B or C>]
     283                          [--negate <yes or no>] [--overwrite] [-w <dir>] -i <infile> -o <outfile>
     284     RDKitFilterPAINS.py -h | --help | -e | --examples
     285
     286 Description:
     287     Filter Pan-assay Interference molecules (PAINS) [ Ref 130 - 131 ] from an input
     288     file by performing a substructure search using SMARTS pattern specified in
     289     MAYACHEMTOOLS/lib/data/PAINSFilter.csv file and write out appropriate
     290     molecules to an output file or simply count the number of filtered molecules.
     291
     292     The supported input file formats are: SD (.sdf, .sd), SMILES (.smi, .csv,
     293      .tsv, .txt)
     294
     295     The supported output file formats are: SD (.sdf, .sd), SMILES (.smi)
     296
     297 Options:
     298     -e, --examples
     299         Print examples.
     300     -h, --help
     301         Print this help message.
     302     -i, --infile <infile>
     303         Input file name.
     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: removeHydrogens,yes,sanitize,yes,strictParsing,yes
     310             SMILES: smilesColumn,1,smilesNameColumn,2,smilesDelimiter,space,
     311                 smilesTitleLine,auto,sanitize,yes
     312
     313         Possible values for smilesDelimiter: space, comma or tab.
     314     -m, --mode <filter or count>  [default: filter]
     315         Specify whether to filter the matched molecules and write out the rest of the
     316         molecules to an outfile or simply count the number of matched molecules
     317         marked for filtering.
     318     -n, --negate <yes or no>  [default: no]
     319         Specify whether to filter molecules not matching the PAINS filters specified by
     320         SMARTS patterns.
     321     -o, --outfile <outfile>
     322         Output file name.
     323     --outfileParams <Name,Value,...>  [default: auto]
     324         A comma delimited list of parameter name and value pairs for writing
     325         molecules to files. The supported parameter names for different file
     326         formats, along with their default values, are shown below:
     327
     328             SD: compute2DCoords,auto,kekulize,no
     329             SMILES: kekulize,no,smilesDelimiter,space, smilesIsomeric,yes,
     330                 smilesTitleLine,yes
     331
     332         Default value for compute2DCoords: yes for SMILES input file; no for all other
     333         file types.
     334     --overwrite
     335         Overwrite existing files.
     336     -p, --painsMode <All, A, B, or C>  [default: All]
     337         Specify PAINS filter family type to used for filtering molecules.
     338     -w, --workingdir <dir>
     339         Location of working directory which defaults to the current directory.
     340
     341 Examples:
     342     To count the number of molecules not containing any substructure corresponding to
     343     PAINS SMARTS patterns and write out a SMILES file, type:
     344
     345         % RDKitFilterPAINS.py -i Sample.smi -o SampleOut.smi
     346
     347     To only count the number of molecules not containing any substructure corresponding
     348     to PAINS SMARTS patterns without writing out any file file, type:
     349
     350         % RDKitFilterPAINS.py -m count -i Sample.sdf -o SampleOut.smi
     351
     352     To count the number of molecules containing any substructure corresponding to
     353     PAINS SMARTS patterns and write out a SD file with computed 2D coordinates,
     354     type:
     355
     356         % RDKitFilterPAINS.py -n yes -i Sample.smi -o SampleOut.sdf
     357
     358     To count the number of molecules not containing any substructure corresponding to
     359     PAINS SMARTS patterns family of Type A in a CSV SMILS file and write out a SD file, type:
     360
     361         % RDKitFilterPAINS.py --painsMode A --infileParams
     362           "smilesDelimiter,comma,smilesTitleLine,yes,smilesColumn,1,
     363           smilesNameColumn,2" --outfileParams "compute2DCoords,yes"
     364           -i SampleSMILES.csv -o SampleOut.sdf
     365
     366 Author:
     367     Manish Sud(msud@san.rr.com)
     368
     369 See also:
     370     RDKitConvertFileFormat.py, RDKitSearchSMARTS.py
     371
     372 Copyright:
     373     Copyright (C) 2018 Manish Sud. All rights reserved.
     374
     375     The functionality available in this script is implemented using RDKit, an
     376     open source toolkit for cheminformatics developed by Greg Landrum.
     377
     378     This file is part of MayaChemTools.
     379
     380     MayaChemTools is free software; you can redistribute it and/or modify it under
     381     the terms of the GNU Lesser General Public License as published by the Free
     382     Software Foundation; either version 3 of the License, or (at your option) any
     383     later version.
     384
     385 """
     386
     387 if __name__ == "__main__":
     388     main()