1
2
3
4
5
6 from __future__ import print_function
7
8 _version = "0.3.2"
9
10 _usage = """
11 ShowFeats [optional args] <filenames>
12
13 if "-" is provided as a filename, data will be read from stdin (the console)
14 """
15
16 _welcomeMessage = "This is ShowFeats version %s" % (_version)
17
18 import math
19
20 from rdkit import RDLogger as logging
21 logger = logging.logger()
22 logger.setLevel(logging.INFO)
23
24 from rdkit import Geometry
25 from rdkit.Chem.Features import FeatDirUtilsRD as FeatDirUtils
26
27 _featColors = {
28 'Donor': (0, 1, 1),
29 'Acceptor': (1, 0, 1),
30 'NegIonizable': (1, 0, 0),
31 'PosIonizable': (0, 0, 1),
32 'ZnBinder': (1, .5, .5),
33 'Aromatic': (1, .8, .2),
34 'LumpedHydrophobe': (.5, .25, 0),
35 'Hydrophobe': (.5, .25, 0),
36 }
37
38
40 if math.fabs(v.x) > tol:
41 res = Geometry.Point3D(v.y, -v.x, 0)
42 elif math.fabs(v.y) > tol:
43 res = Geometry.Point3D(-v.y, v.x, 0)
44 elif math.fabs(v.z) > tol:
45 res = Geometry.Point3D(1, 0, 0)
46 else:
47 raise ValueError('cannot find normal to the null vector')
48 res.Normalize()
49 return res
50
51
52 _canonArrowhead = None
53
54
56 global _canonArrowhead
57 startP = RDGeometry.Point3D(0, 0, headFrac)
58 _canonArrowhead = [startP]
59
60 scale = headFrac * aspect
61 baseV = RDGeometry.Point3D(scale, 0, 0)
62 _canonArrowhead.append(baseV)
63
64 twopi = 2 * math.pi
65 for i in range(1, nSteps):
66 v = RDGeometry.Point3D(scale * math.cos(i * twopi), scale * math.sin(i * twopi), 0)
67 _canonArrowhead.append(v)
68
69
70 _globalArrowCGO = []
71 _globalSphereCGO = []
72
73 BEGIN = 2
74 END = 3
75 TRIANGLE_FAN = 6
76 COLOR = 6
77 VERTEX = 4
78 NORMAL = 5
79 SPHERE = 7
80 CYLINDER = 9
81 ALPHA = 25
82
83
84 -def _cgoArrowhead(viewer, tail, head, radius, color, label, headFrac=0.3, nSteps=10, aspect=.5):
85 global _globalArrowCGO
86 delta = head - tail
87 normal = _getVectNormal(delta)
88 delta.Normalize()
89
90 dv = head - tail
91 dv.Normalize()
92 dv *= headFrac
93 startP = head
94
95 normal *= headFrac * aspect
96
97 cgo = [BEGIN, TRIANGLE_FAN, COLOR, color[0], color[1], color[2], NORMAL, dv.x, dv.y, dv.z, VERTEX,
98 head.x + dv.x, head.y + dv.y, head.z + dv.z]
99 base = [BEGIN, TRIANGLE_FAN, COLOR, color[0], color[1], color[2], NORMAL, -dv.x, -dv.y, -dv.z,
100 VERTEX, head.x, head.y, head.z]
101 v = startP + normal
102 cgo.extend([NORMAL, normal.x, normal.y, normal.z])
103 cgo.extend([VERTEX, v.x, v.y, v.z])
104 base.extend([VERTEX, v.x, v.y, v.z])
105 for i in range(1, nSteps):
106 v = FeatDirUtils.ArbAxisRotation(360. / nSteps * i, delta, normal)
107 cgo.extend([NORMAL, v.x, v.y, v.z])
108 v += startP
109 cgo.extend([VERTEX, v.x, v.y, v.z])
110 base.extend([VERTEX, v.x, v.y, v.z])
111
112 cgo.extend([NORMAL, normal.x, normal.y, normal.z])
113 cgo.extend([VERTEX, startP.x + normal.x, startP.y + normal.y, startP.z + normal.z])
114 base.extend([VERTEX, startP.x + normal.x, startP.y + normal.y, startP.z + normal.z])
115 cgo.append(END)
116 base.append(END)
117 cgo.extend(base)
118
119
120 _globalArrowCGO.extend(cgo)
121
122
123 -def ShowArrow(viewer, tail, head, radius, color, label, transparency=0, includeArrowhead=True):
124 global _globalArrowCGO
125 if transparency:
126 _globalArrowCGO.extend([ALPHA, 1 - transparency])
127 else:
128 _globalArrowCGO.extend([ALPHA, 1])
129 _globalArrowCGO.extend([CYLINDER,
130 tail.x,
131 tail.y,
132 tail.z,
133 head.x,
134 head.y,
135 head.z,
136 radius * .10,
137 color[0],
138 color[1],
139 color[2],
140 color[0],
141 color[1],
142 color[2], ])
143
144 if includeArrowhead:
145 _cgoArrowhead(viewer, tail, head, radius, color, label)
146
147
148 -def ShowMolFeats(mol, factory, viewer, radius=0.5, confId=-1, showOnly=True, name='',
149 transparency=0.0, colors=None, excludeTypes=[], useFeatDirs=True, featLabel=None,
150 dirLabel=None, includeArrowheads=True, writeFeats=False, showMol=True,
151 featMapFile=False):
152 global _globalSphereCGO
153 if not name:
154 if mol.HasProp('_Name'):
155 name = mol.GetProp('_Name')
156 else:
157 name = 'molecule'
158 if not colors:
159 colors = _featColors
160
161 if showMol:
162 viewer.ShowMol(mol, name=name, showOnly=showOnly, confId=confId)
163
164 molFeats = factory.GetFeaturesForMol(mol)
165 if not featLabel:
166 featLabel = '%s-feats' % name
167 viewer.server.resetCGO(featLabel)
168 if not dirLabel:
169 dirLabel = featLabel + "-dirs"
170 viewer.server.resetCGO(dirLabel)
171
172 for i, feat in enumerate(molFeats):
173 family = feat.GetFamily()
174 if family in excludeTypes:
175 continue
176 pos = feat.GetPos(confId)
177 color = colors.get(family, (.5, .5, .5))
178 nm = '%s(%d)' % (family, i + 1)
179
180 if transparency:
181 _globalSphereCGO.extend([ALPHA, 1 - transparency])
182 else:
183 _globalSphereCGO.extend([ALPHA, 1])
184 _globalSphereCGO.extend([COLOR, color[0], color[1], color[2], SPHERE, pos.x, pos.y, pos.z,
185 radius])
186 if writeFeats:
187 aidText = ' '.join([str(x + 1) for x in feat.GetAtomIds()])
188 print('%s\t%.3f\t%.3f\t%.3f\t1.0\t# %s' % (family, pos.x, pos.y, pos.z, aidText))
189
190 if featMapFile:
191 print(" family=%s pos=(%.3f,%.3f,%.3f) weight=1.0" % (family, pos.x, pos.y, pos.z), end='',
192 file=featMapFile)
193
194 if useFeatDirs:
195 ps = []
196 if family == 'Aromatic':
197 ps, fType = FeatDirUtils.GetAromaticFeatVects(
198 mol.GetConformer(confId), feat.GetAtomIds(), pos, scale=1.0)
199 elif family == 'Donor':
200 aids = feat.GetAtomIds()
201 if len(aids) == 1:
202 featAtom = mol.GetAtomWithIdx(aids[0])
203 hvyNbrs = [x for x in featAtom.GetNeighbors() if x.GetAtomicNum() != 1]
204 if len(hvyNbrs) == 1:
205 ps, fType = FeatDirUtils.GetDonor1FeatVects(mol.GetConformer(confId), aids, scale=1.0)
206 elif len(hvyNbrs) == 2:
207 ps, fType = FeatDirUtils.GetDonor2FeatVects(mol.GetConformer(confId), aids, scale=1.0)
208 elif len(hvyNbrs) == 3:
209 ps, fType = FeatDirUtils.GetDonor3FeatVects(mol.GetConformer(confId), aids, scale=1.0)
210 elif family == 'Acceptor':
211 aids = feat.GetAtomIds()
212 if len(aids) == 1:
213 featAtom = mol.GetAtomWithIdx(aids[0])
214 hvyNbrs = [x for x in featAtom.GetNeighbors() if x.GetAtomicNum() != 1]
215 if len(hvyNbrs) == 1:
216 ps, fType = FeatDirUtils.GetAcceptor1FeatVects(
217 mol.GetConformer(confId), aids, scale=1.0)
218 elif len(hvyNbrs) == 2:
219 ps, fType = FeatDirUtils.GetAcceptor2FeatVects(
220 mol.GetConformer(confId), aids, scale=1.0)
221 elif len(hvyNbrs) == 3:
222 ps, fType = FeatDirUtils.GetAcceptor3FeatVects(
223 mol.GetConformer(confId), aids, scale=1.0)
224
225 for tail, head in ps:
226 ShowArrow(viewer, tail, head, radius, color, dirLabel, transparency=transparency,
227 includeArrowhead=includeArrowheads)
228 if featMapFile:
229 vect = head - tail
230 print('dir=(%.3f,%.3f,%.3f)' % (vect.x, vect.y, vect.z), end='', file=featMapFile)
231
232 if featMapFile:
233 aidText = ' '.join([str(x + 1) for x in feat.GetAtomIds()])
234 print('# %s' % (aidText), file=featMapFile)
235
236
237 import sys, os, getopt
238 from rdkit import RDConfig
239 from optparse import OptionParser
240 parser = OptionParser(_usage, version='%prog ' + _version)
241
242 parser.add_option('-x', '--exclude', default='',
243 help='provide a list of feature names that should be excluded')
244 parser.add_option('-f', '--fdef', default=os.path.join(RDConfig.RDDataDir, 'BaseFeatures.fdef'),
245 help='provide the name of the feature definition (fdef) file.')
246 parser.add_option('--noDirs', '--nodirs', dest='useDirs', default=True, action='store_false',
247 help='do not draw feature direction indicators')
248 parser.add_option('--noHeads', dest='includeArrowheads', default=True, action='store_false',
249 help='do not draw arrowheads on the feature direction indicators')
250 parser.add_option('--noClear', '--noClear', dest='clearAll', default=False, action='store_true',
251 help='do not clear PyMol on startup')
252 parser.add_option('--noMols', '--nomols', default=False, action='store_true',
253 help='do not draw the molecules')
254 parser.add_option('--writeFeats', '--write', default=False, action='store_true',
255 help='print the feature information to the console')
256 parser.add_option('--featMapFile', '--mapFile', default='',
257 help='save a feature map definition to the specified file')
258 parser.add_option('--verbose', default=False, action='store_true', help='be verbose')
259
260 if __name__ == '__main__':
261 from rdkit import Chem
262 from rdkit.Chem import AllChem
263 from rdkit.Chem.PyMol import MolViewer
264
265 options, args = parser.parse_args()
266 if len(args) < 1:
267 parser.error('please provide either at least one sd or mol file')
268
269 try:
270 v = MolViewer()
271 except Exception:
272 logger.error(
273 'Unable to connect to PyMol server.\nPlease run ~landrgr1/extern/PyMol/launch.sh to start it.')
274 sys.exit(1)
275 if options.clearAll:
276 v.DeleteAll()
277
278 try:
279 fdef = open(options.fdef, 'r').read()
280 except IOError:
281 logger.error('ERROR: Could not open fdef file %s' % options.fdef)
282 sys.exit(1)
283
284 factory = AllChem.BuildFeatureFactoryFromString(fdef)
285
286 if options.writeFeats:
287 print('# Family \tX \tY \tZ \tRadius\t # Atom_ids')
288
289 if options.featMapFile:
290 if options.featMapFile == '-':
291 options.featMapFile = sys.stdout
292 else:
293 options.featMapFile = file(options.featMapFile, 'w+')
294 print('# Feature map generated by ShowFeats v%s' % _version, file=options.featMapFile)
295 print("ScoreMode=All", file=options.featMapFile)
296 print("DirScoreMode=Ignore", file=options.featMapFile)
297 print("BeginParams", file=options.featMapFile)
298 for family in factory.GetFeatureFamilies():
299 print(" family=%s width=1.0 radius=3.0" % family, file=options.featMapFile)
300 print("EndParams", file=options.featMapFile)
301 print("BeginPoints", file=options.featMapFile)
302
303 i = 1
304 for midx, molN in enumerate(args):
305 if molN != '-':
306 featLabel = '%s_Feats' % molN
307 else:
308 featLabel = 'Mol%d_Feats' % (midx + 1)
309
310 v.server.resetCGO(featLabel)
311
312 v.server.sphere((0, 0, 0), .01, (1, 0, 1), featLabel)
313 dirLabel = featLabel + "-dirs"
314
315 v.server.resetCGO(dirLabel)
316
317 v.server.cylinder((0, 0, 0), (.01, .01, .01), .01, (1, 0, 1), dirLabel)
318
319 if molN != '-':
320 try:
321 ms = Chem.SDMolSupplier(molN)
322 except Exception:
323 logger.error('Problems reading input file: %s' % molN)
324 ms = []
325 else:
326 ms = Chem.SDMolSupplier()
327 ms.SetData(sys.stdin.read())
328
329 for m in ms:
330 nm = 'Mol_%d' % (i)
331 if m.HasProp('_Name'):
332 nm += '_' + m.GetProp('_Name')
333 if options.verbose:
334 if m.HasProp('_Name'):
335 print("#Molecule: %s" % m.GetProp('_Name'))
336 else:
337 print("#Molecule: %s" % nm)
338 ShowMolFeats(m, factory, v, transparency=0.25, excludeTypes=options.exclude, name=nm,
339 showOnly=False, useFeatDirs=options.useDirs, featLabel=featLabel,
340 dirLabel=dirLabel, includeArrowheads=options.includeArrowheads,
341 writeFeats=options.writeFeats, showMol=not options.noMols,
342 featMapFile=options.featMapFile)
343 i += 1
344 if not i % 100:
345 logger.info("Done %d poses" % i)
346 if ms:
347 v.server.renderCGO(_globalSphereCGO, featLabel, 1)
348 if options.useDirs:
349 v.server.renderCGO(_globalArrowCGO, dirLabel, 1)
350
351 if options.featMapFile:
352 print("EndPoints", file=options.featMapFile)
353
354 sys.exit(0)
355