1
2
3
4
5
6 from __future__ import print_function
7
8 _version = "0.3.2"
9
10
11 _usage="""
12 ShowFeats [optional args] <filenames>
13
14 if "-" is provided as a filename, data will be read from stdin (the console)
15 """
16
17 _welcomeMessage="This is ShowFeats version %s"%(_version)
18
19
20 import math
21
22 from rdkit import RDLogger as logging
23 logger = logging.logger()
24 logger.setLevel(logging.INFO)
25
26 from rdkit import Geometry
27 from rdkit.Chem.Features import FeatDirUtilsRD as FeatDirUtils
28
29 _featColors = {
30 'Donor':(0,1,1),
31 'Acceptor':(1,0,1),
32 'NegIonizable':(1,0,0),
33 'PosIonizable':(0,0,1),
34 'ZnBinder':(1,.5,.5),
35 'Aromatic':(1,.8,.2),
36 'LumpedHydrophobe':(.5,.25,0),
37 'Hydrophobe':(.5,.25,0),
38 }
39
41 if math.fabs(v.x)>tol:
42 res = Geometry.Point3D(v.y,-v.x,0)
43 elif math.fabs(v.y)>tol:
44 res = Geometry.Point3D(-v.y,v.x,0)
45 elif math.fabs(v.z)>tol:
46 res = Geometry.Point3D(1,0,0)
47 else:
48 raise ValueError('cannot find normal to the null vector')
49 res.Normalize()
50 return res
51
52 _canonArrowhead=None
54 global _canonArrowhead
55 startP = RDGeometry.Point3D(0,0,headFrac)
56 _canonArrowhead=[startP]
57
58 scale = headFrac*aspect
59 baseV = RDGeometry.Point3D(scale,0,0)
60 _canonArrowhead.append(baseV)
61
62 twopi = 2*math.pi
63 for i in range(1,nSteps):
64 v = RDGeometry.Point3D(scale*math.cos(i*twopi),scale*math.sin(i*twopi),0)
65 _canonArrowhead.append(v)
66
67
68 _globalArrowCGO=[]
69 _globalSphereCGO=[]
70
71 BEGIN=2
72 END=3
73 TRIANGLE_FAN=6
74 COLOR=6
75 VERTEX=4
76 NORMAL=5
77 SPHERE=7
78 CYLINDER=9
79 ALPHA=25
80
81 -def _cgoArrowhead(viewer,tail,head,radius,color,label,headFrac=0.3,nSteps=10,aspect=.5):
82 global _globalArrowCGO
83 delta = head-tail
84 normal = _getVectNormal(delta)
85 delta.Normalize()
86
87 dv = head-tail
88 dv.Normalize()
89 dv *= headFrac
90 startP = head
91
92 normal*=headFrac*aspect
93
94 cgo = [BEGIN,TRIANGLE_FAN,
95 COLOR,color[0],color[1],color[2],
96 NORMAL,dv.x,dv.y,dv.z,
97 VERTEX,head.x+dv.x,head.y+dv.y,head.z+dv.z]
98 base = [BEGIN,TRIANGLE_FAN,
99 COLOR,color[0],color[1],color[2],
100 NORMAL,-dv.x,-dv.y,-dv.z,
101 VERTEX,head.x,head.y,head.z]
102 v = startP+normal
103 cgo.extend([NORMAL,normal.x,normal.y,normal.z])
104 cgo.extend([VERTEX,v.x,v.y,v.z])
105 base.extend([VERTEX,v.x,v.y,v.z])
106 for i in range(1,nSteps):
107 v = FeatDirUtils.ArbAxisRotation(360./nSteps*i,delta,normal)
108 cgo.extend([NORMAL,v.x,v.y,v.z])
109 v += startP
110 cgo.extend([VERTEX,v.x,v.y,v.z])
111 base.extend([VERTEX,v.x,v.y,v.z])
112
113 cgo.extend([NORMAL,normal.x,normal.y,normal.z])
114 cgo.extend([VERTEX,startP.x+normal.x,startP.y+normal.y,startP.z+normal.z])
115 base.extend([VERTEX,startP.x+normal.x,startP.y+normal.y,startP.z+normal.z])
116 cgo.append(END)
117 base.append(END)
118 cgo.extend(base)
119
120
121 _globalArrowCGO.extend(cgo)
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,tail.x,tail.y,tail.z,
130 head.x,head.y,head.z,
131 radius*.10,
132 color[0],color[1],color[2],
133 color[0],color[1],color[2],
134 ])
135
136 if includeArrowhead:
137 _cgoArrowhead(viewer,tail,head,radius,color,label)
138
139
140 -def ShowMolFeats(mol,factory,viewer,radius=0.5,confId=-1,showOnly=True,
141 name='',transparency=0.0,colors=None,excludeTypes=[],
142 useFeatDirs=True,featLabel=None,dirLabel=None,includeArrowheads=True,
143 writeFeats=False,showMol=True,featMapFile=False):
144 global _globalSphereCGO
145 if not name:
146 if mol.HasProp('_Name'):
147 name = mol.GetProp('_Name')
148 else:
149 name = 'molecule'
150 if not colors:
151 colors = _featColors
152
153 if showMol:
154 viewer.ShowMol(mol,name=name,showOnly=showOnly,confId=confId)
155
156 molFeats=factory.GetFeaturesForMol(mol)
157 if not featLabel:
158 featLabel='%s-feats'%name
159 viewer.server.resetCGO(featLabel)
160 if not dirLabel:
161 dirLabel=featLabel+"-dirs"
162 viewer.server.resetCGO(dirLabel)
163
164 for i,feat in enumerate(molFeats):
165 family=feat.GetFamily()
166 if family in excludeTypes:
167 continue
168 pos = feat.GetPos(confId)
169 color = colors.get(family,(.5,.5,.5))
170 nm = '%s(%d)'%(family,i+1)
171
172 if transparency:
173 _globalSphereCGO.extend([ALPHA,1-transparency])
174 else:
175 _globalSphereCGO.extend([ALPHA,1])
176 _globalSphereCGO.extend([COLOR,color[0],color[1],color[2],
177 SPHERE,pos.x,pos.y,pos.z,
178 radius])
179 if writeFeats:
180 aidText = ' '.join([str(x+1) for x in feat.GetAtomIds()])
181 print('%s\t%.3f\t%.3f\t%.3f\t1.0\t# %s'%(family,pos.x,pos.y,pos.z,aidText))
182
183 if featMapFile:
184 print(" family=%s pos=(%.3f,%.3f,%.3f) weight=1.0"%(family,pos.x,pos.y,pos.z),end='',file=featMapFile)
185
186 if useFeatDirs:
187 ps = []
188 if family=='Aromatic':
189 ps,fType = FeatDirUtils.GetAromaticFeatVects(mol.GetConformer(confId),
190 feat.GetAtomIds(),pos,
191 scale=1.0)
192 elif family=='Donor':
193 aids = feat.GetAtomIds()
194 if len(aids)==1:
195 featAtom=mol.GetAtomWithIdx(aids[0])
196 hvyNbrs=[x for x in featAtom.GetNeighbors() if x.GetAtomicNum()!=1]
197 if len(hvyNbrs)==1:
198 ps,fType = FeatDirUtils.GetDonor1FeatVects(mol.GetConformer(confId),
199 aids,scale=1.0)
200 elif len(hvyNbrs)==2:
201 ps,fType = FeatDirUtils.GetDonor2FeatVects(mol.GetConformer(confId),
202 aids,scale=1.0)
203 elif len(hvyNbrs)==3:
204 ps,fType = FeatDirUtils.GetDonor3FeatVects(mol.GetConformer(confId),
205 aids,scale=1.0)
206 elif family=='Acceptor':
207 aids = feat.GetAtomIds()
208 if len(aids)==1:
209 featAtom=mol.GetAtomWithIdx(aids[0])
210 hvyNbrs=[x for x in featAtom.GetNeighbors() if x.GetAtomicNum()!=1]
211 if len(hvyNbrs)==1:
212 ps,fType = FeatDirUtils.GetAcceptor1FeatVects(mol.GetConformer(confId),
213 aids,scale=1.0)
214 elif len(hvyNbrs)==2:
215 ps,fType = FeatDirUtils.GetAcceptor2FeatVects(mol.GetConformer(confId),
216 aids,scale=1.0)
217 elif len(hvyNbrs)==3:
218 ps,fType = FeatDirUtils.GetAcceptor3FeatVects(mol.GetConformer(confId),
219 aids,scale=1.0)
220
221 for tail,head in ps:
222 ShowArrow(viewer,tail,head,radius,color,dirLabel,
223 transparency=transparency,includeArrowhead=includeArrowheads)
224 if featMapFile:
225 vect = head-tail
226 print('dir=(%.3f,%.3f,%.3f)'%(vect.x,vect.y,vect.z),end='',file=featMapFile)
227
228 if featMapFile:
229 aidText = ' '.join([str(x+1) for x in feat.GetAtomIds()])
230 print('# %s'%(aidText),file=featMapFile)
231
232
233
234
235
236 import sys,os,getopt
237 from rdkit import RDConfig
238 from optparse import OptionParser
239 parser=OptionParser(_usage,version='%prog '+_version)
240
241 parser.add_option('-x','--exclude',default='',
242 help='provide a list of feature names that should be excluded')
243 parser.add_option('-f','--fdef',default=os.path.join(RDConfig.RDDataDir,'BaseFeatures.fdef'),
244 help='provide the name of the feature definition (fdef) file.')
245 parser.add_option('--noDirs','--nodirs',dest='useDirs',default=True,action='store_false',
246 help='do not draw feature direction indicators')
247 parser.add_option('--noHeads',dest='includeArrowheads',default=True,action='store_false',
248 help='do not draw arrowheads on the feature direction indicators')
249 parser.add_option('--noClear','--noClear',dest='clearAll',default=False,action='store_true',
250 help='do not clear PyMol on startup')
251 parser.add_option('--noMols','--nomols',default=False,action='store_true',
252 help='do not draw the molecules')
253 parser.add_option('--writeFeats','--write',default=False,action='store_true',
254 help='print the feature information to the console')
255 parser.add_option('--featMapFile','--mapFile',default='',
256 help='save a feature map definition to the specified file')
257 parser.add_option('--verbose',default=False,action='store_true',
258 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('Unable to connect to PyMol server.\nPlease run ~landrgr1/extern/PyMol/launch.sh to start it.')
273 sys.exit(1)
274 if options.clearAll:
275 v.DeleteAll()
276
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,
340 useFeatDirs=options.useDirs,
341 featLabel=featLabel,dirLabel=dirLabel,
342 includeArrowheads=options.includeArrowheads,
343 writeFeats=options.writeFeats,showMol=not options.noMols,
344 featMapFile=options.featMapFile)
345 i += 1
346 if not i%100:
347 logger.info("Done %d poses"%i)
348 if ms:
349 v.server.renderCGO(_globalSphereCGO,featLabel,1)
350 if options.useDirs:
351 v.server.renderCGO(_globalArrowCGO,dirLabel,1)
352
353 if options.featMapFile:
354 print("EndPoints",file=options.featMapFile)
355
356
357 sys.exit(0)
358