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