Package rdkit :: Package Chem :: Package Features :: Module ShowFeats
[hide private]
[frames] | no frames]

Source Code for Module rdkit.Chem.Features.ShowFeats

  1  # $Id: ShowFeats.py 537 2007-08-20 14:54:35Z landrgr1 $ 
  2  # 
  3  # Created by Greg Landrum Aug 2006 
  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  #set up the logger: 
 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   
40 -def _getVectNormal(v,tol=1e-4):
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
53 -def _buildCanonArrowhead(headFrac,nSteps,aspect):
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 # taken from pymol's cgo.py 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 #viewer.server.renderCGO(cgo,label) 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 # this is a big of kludgery to work around what seems to be a pymol cgo bug: 312 v.server.sphere((0,0,0),.01,(1,0,1),featLabel) 313 dirLabel=featLabel+"-dirs" 314 315 v.server.resetCGO(dirLabel) 316 # this is a big of kludgery to work around what seems to be a pymol cgo bug: 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