1#
2#  Copyright (C) 2005-2006 Rational Discovery LLC
3#
4#   @@ All Rights Reserved @@
5#  This file is part of the RDKit.
6#  The contents are covered by the terms of the BSD license
7#  which is included in the file license.txt, found at the root
8#  of the RDKit source tree.
9#
10
11
12import argparse
13import re
14import os
15
16from rdkit import Chem
17from rdkit import RDLogger
18from rdkit.Chem import ChemicalFeatures
19
20logger = RDLogger.logger()
21splitExpr = re.compile(r'[ \t,]')
22
23
24def GetAtomFeatInfo(factory, mol):
25  res = [None] * mol.GetNumAtoms()
26  feats = factory.GetFeaturesForMol(mol)
27  for feat in feats:
28    ids = feat.GetAtomIds()
29    feature = "%s-%s" % (feat.GetFamily(), feat.GetType())
30    for id_ in ids:
31      if res[id_] is None:
32        res[id_] = []
33      res[id_].append(feature)
34  return res
35
36
37def initParser():
38  """ Initialize the parser """
39  parser = argparse.ArgumentParser(description='Determine pharmacophore features of molecules',
40                                   epilog=_splashMessage,
41                                   formatter_class=argparse.RawDescriptionHelpFormatter)
42  parser.add_argument('-r', dest='reverseIt', default=False, action='store_true',
43                      help='Set to get atoms lists for each feature.')
44  parser.add_argument('-n', dest='maxLines', default=-1, help=argparse.SUPPRESS, type=int)
45  parser.add_argument('fdefFilename', type=existingFile,
46                      help='Pharmacophore feature definition file')
47  parser.add_argument('smilesFilename', type=existingFile,
48                      help='The smiles file should have SMILES in the first column')
49  return parser
50
51
52_splashMessage = """
53-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
54  FeatFinderCLI
55  Part of the RDKit (http://www.rdkit.org)
56-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
57"""
58
59
60def existingFile(filename):
61  """ 'type' for argparse - check that filename exists """
62  if not os.path.exists(filename):
63    raise argparse.ArgumentTypeError("{0} does not exist".format(filename))
64  return filename
65
66
67def processArgs(args, parser):
68  try:
69    factory = ChemicalFeatures.BuildFeatureFactory(args.fdefFilename)
70  except Exception:
71    parser.error("Could not parse Fdef file {0.fdefFilename}.".format(args))
72
73  with open(args.smilesFilename) as inF:
74    for lineNo, line in enumerate(inF, 1):
75      if lineNo == args.maxLines + 1:
76        break
77      smi = splitExpr.split(line.strip())[0].strip()
78      mol = Chem.MolFromSmiles(smi)
79      if mol is None:
80        logger.warning("Could not process smiles '%s' on line %d." % (smi, lineNo))
81        continue
82
83      print('Mol-%d\t%s' % (lineNo, smi))
84      if args.reverseIt:
85        feats = factory.GetFeaturesForMol(mol)
86        for feat in feats:
87          print('\t%s-%s: ' % (feat.GetFamily(), feat.GetType()), end='')
88          print(', '.join([str(x) for x in feat.GetAtomIds()]))
89      else:
90        featInfo = GetAtomFeatInfo(factory, mol)
91        for i, v in enumerate(featInfo):
92          print('\t% 2s(%d)' % (mol.GetAtomWithIdx(i).GetSymbol(), i + 1), end='')
93          if v:
94            print('\t', ', '.join(v))
95          else:
96            print()
97
98
99def main():
100  """ Main application """
101  parser = initParser()
102  args = parser.parse_args()
103  processArgs(args, parser)
104
105
106if __name__ == '__main__':
107  main()
108