1
2
3
4
5
6
7
8
9
10 from __future__ import print_function
11
12 import argparse
13 import re
14 import os
15
16 from rdkit import Chem
17 from rdkit import RDLogger
18 from rdkit.Chem import ChemicalFeatures
19
20 logger = RDLogger.logger()
21 splitExpr = re.compile(r'[ \t,]')
22
23
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
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
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
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
100 """ Main application """
101 parser = initParser()
102 args = parser.parse_args()
103 processArgs(args, parser)
104
105
106 if __name__ == '__main__':
107 main()
108