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