Source Code for Module Chem.Subshape.SubshapeBuilder

  1  # $Id: SubshapeBuilder.py 746 2008-07-07 13:21:24Z glandrum $ 
  2  # 
  3  # Copyright (C) 2007 by Greg Landrum  
  4  #  All rights reserved 
  5  # 
  6  import Chem,Geometry 
  7  from Chem import AllChem 
  8  from Chem.Subshape import SubshapeObjects 
  9  from Chem.Subshape import BuilderUtils 
 10  import time,cPickle 
 11   
 12  #----------------------------------------------------------------------------- 
 13 -class SubshapeCombineOperations(object): 
 14    UNION=0 
 15    SUM=1 
 16    INTERSECT=2 
 17   
 18  #----------------------------------------------------------------------------- 
 19 -class SubshapeBuilder(object): 
 20    gridDims=(20,15,10) 
 21    gridSpacing=0.5 
 22    winRad=3.0 
 23    nbrCount=7 
 24    terminalPtRadScale=0.75 
 25    fraction=0.25 
 26    stepSize=1.0 
 27    featFactory=None 
 28   
 29 -  def SampleSubshape(self,subshape1,newSpacing): 
 30      ogrid=subshape1.grid 
 31      rgrid = Geometry.UniformGrid3D(self.gridDims[0],self.gridDims[1],self.gridDims[2], 
 32                                        newSpacing) 
 33      for idx in range(rgrid.GetSize()): 
 34        l = rgrid.GetGridPointLoc(idx) 
 35        v = ogrid.GetValPoint(l) 
 36        rgrid.SetVal(idx,v) 
 37         
 38      res = SubshapeObjects.ShapeWithSkeleton() 
 39      res.grid = rgrid 
 40      return res; 
 41   
 42 -  def GenerateSubshapeShape(self,cmpd,confId=-1,addSkeleton=True,**kwargs): 
 43      shape = SubshapeObjects.ShapeWithSkeleton() 
 44      shape.grid=Geometry.UniformGrid3D(self.gridDims[0],self.gridDims[1],self.gridDims[2], 
 45                                        self.gridSpacing) 
 46      AllChem.EncodeShape(cmpd,shape.grid,ignoreHs=False,confId=confId) 
 47      if addSkeleton: 
 48        conf = cmpd.GetConformer(confId) 
 49        self.GenerateSubshapeSkeleton(shape,conf,kwargs) 
 50      return shape 
 51 -  def __call__(self,cmpd,**kwargs): 
 52      return self.GenerateSubshapeShape(cmpd,**kwargs) 
 53     
 54 -  def GenerateSubshapeSkeleton(self,shape,conf=None,terminalPtsOnly=False,skelFromConf=True): 
 55      if conf and skelFromConf: 
 56        pts = BuilderUtils.FindTerminalPtsFromConformer(conf,self.winRad,self.nbrCount) 
 57      else: 
 58        pts = BuilderUtils.FindTerminalPtsFromShape(shape,self.winRad,self.fraction) 
 59      pts = BuilderUtils.ClusterTerminalPts(pts,self.winRad,self.terminalPtRadScale) 
 60      BuilderUtils.ExpandTerminalPts(shape,pts,self.winRad) 
 61      if len(pts)<3: 
 62        raise ValueError,'only found %d terminals, need at least 3'%len(pts) 
 63       
 64      if not terminalPtsOnly: 
 65        pts = BuilderUtils.AppendSkeletonPoints(shape.grid,pts,self.winRad,self.stepSize) 
 66      for i,pt in enumerate(pts): 
 67        BuilderUtils.CalculateDirectionsAtPoint(pt,shape.grid,self.winRad) 
 68      if conf and self.featFactory: 
 69        BuilderUtils.AssignMolFeatsToPoints(pts,conf.GetOwningMol(),self.featFactory,self.winRad) 
 70      shape.skelPts=pts 
 71   
 72 -  def CombineSubshapes(self,subshape1,subshape2,operation=SubshapeCombineOperations.UNION): 
 73      import copy 
 74      cs = copy.deepcopy(subshape1) 
 75      if operation==SubshapeCombineOperations.UNION: 
 76        cs.grid |= subshape2.grid 
 77      elif operation==SubshapeCombineOperations.SUM: 
 78        cs.grid +=  subshape2.grid 
 79      elif operation==SubshapeCombineOperations.INTERSECT: 
 80        cs.grid &=  subshape2.grid 
 81      else: 
 82        raise ValueError,'bad combination operation' 
 83      return cs 
 84   
 85       
 86  if __name__=='__main__': 
 87    from Chem import AllChem,ChemicalFeatures 
 88    from Chem.PyMol import MolViewer 
 89    #cmpd = Chem.MolFromSmiles('CCCc1cc(C(=O)O)ccc1') 
 90    #cmpd = Chem.AddHs(cmpd) 
 91    if 1: 
 92      cmpd = Chem.MolFromSmiles('C1=CC=C1C#CC1=CC=C1') 
 93      cmpd = Chem.AddHs(cmpd) 
 94      AllChem.EmbedMolecule(cmpd) 
 95      AllChem.UFFOptimizeMolecule(cmpd) 
 96      AllChem.CanonicalizeMol(cmpd) 
 97      print >>file('testmol.mol','w+'),Chem.MolToMolBlock(cmpd) 
 98    else: 
 99      cmpd = Chem.MolFromMolFile('testmol.mol') 
100    builder=SubshapeBuilder() 
101    if 1: 
102      shape=builder.GenerateSubshapeShape(cmpd) 
103    v = MolViewer() 
104    if 1: 
105      import tempfile 
106      tmpFile = tempfile.mktemp('.grd') 
107      v.server.deleteAll() 
108      Geometry.WriteGridToFile(shape.grid,tmpFile) 
109      time.sleep(1) 
110      v.ShowMol(cmpd,name='testMol',showOnly=True) 
111      v.server.loadSurface(tmpFile,'testGrid','',2.5) 
112    v.server.resetCGO('*') 
113   
114    cPickle.dump(shape,file('subshape.pkl','w+')) 
115    for i,pt in enumerate(shape.skelPts): 
116      v.server.sphere(tuple(pt.location),.5,(1,0,1),'Pt-%d'%i) 
117      if not hasattr(pt,'shapeDirs'): continue 
118      momBeg = pt.location-pt.shapeDirs[0] 
119      momEnd = pt.location+pt.shapeDirs[0] 
120      v.server.cylinder(tuple(momBeg),tuple(momEnd),.1,(1,0,1),'v-%d'%i) 
121