Source Code for Module Chem.Pharm3D.Pharmacophore

  1  # $Id: Pharmacophore.py 746 2008-07-07 13:21:24Z glandrum $ 
  2  # 
  3  # Copyright (C) 2004-2008 Greg Landrum and Rational Discovery LLC 
  4  # 
  5  #   @@ All Rights Reserved  @@ 
  6  # 
  7  import Geometry 
  8  import numpy 
  9  import Chem 
 10  from Chem import ChemicalFeatures 
 11   
 12 -class Pharmacophore: 
 13 -  def __init__(self, feats,initMats=True): 
 14      self._initializeFeats(feats) 
 15      nf = len(feats) 
 16      self._boundsMat = numpy.zeros((nf, nf), numpy.float) 
 17      self._boundsMat2D = numpy.zeros((nf, nf), numpy.int) 
 18      if initMats: 
 19        self._initializeMatrices() 
 20   
 21 -  def _initializeFeats(self,feats): 
 22      self._feats = [] 
 23      for feat in feats: 
 24        if isinstance(feat,ChemicalFeatures.MolChemicalFeature): 
 25          pos = feat.GetPos() 
 26          newFeat = ChemicalFeatures.FreeChemicalFeature(feat.GetFamily(),feat.GetType(), 
 27                                             Geometry.Point3D(pos[0],pos[1],pos[2])) 
 28          self._feats.append(newFeat) 
 29        else: 
 30          self._feats.append(feat) 
 31   
 32 -  def _initializeMatrices(self): 
 33      # initialize the bounds matrix with distances to start with 
 34      nf = len(self._feats) 
 35      for i in range(1, nf): 
 36        loci = self._feats[i].GetPos() 
 37        for j in range(i): 
 38          locj = self._feats[j].GetPos() 
 39          dist = loci.Distance(locj) 
 40          self._boundsMat[i,j] = dist 
 41          self._boundsMat[j,i] = dist 
 42      for i in range(nf): 
 43        for j in range(i+1,nf): 
 44          self._boundsMat2D[i,j] = 1000 
 45   
 46   
 47 -  def getFeatures(self): 
 48      return self._feats 
 49   
 50 -  def getFeature(self, i): 
 51      return self._feats[i] 
 52   
 53 -  def getUpperBound(self, i, j): 
 54      if (i > j): 
 55        j,i = i,j 
 56      return self._boundsMat[i, j] 
 57   
 58 -  def getLowerBound(self, i, j): 
 59      if (j > i): 
 60        j,i = i,j 
 61      return self._boundsMat[i,j] 
 62   
 63 -  def _checkBounds(self,i,j): 
 64      " raises ValueError on failure " 
 65      nf = len(self._feats) 
 66      if (i < 0) or (i >= nf): 
 67        raise ValueError, "Index out of bound" 
 68      if (j < 0) or (j >= nf): 
 69        raise ValueError, "Index out of bound" 
 70      return True 
 71   
 72 -  def setUpperBound(self, i, j, val, checkBounds=False): 
 73      if (checkBounds): self._checkBounds(i,j) 
 74      if (i > j): 
 75        j,i = i,j 
 76      self._boundsMat[i,j] = val 
 77   
 78 -  def setLowerBound(self, i, j, val, checkBounds=False): 
 79      if (checkBounds): self._checkBounds(i,j) 
 80      if (j > i): 
 81        j,i = i,j 
 82      self._boundsMat[i,j] = val 
 83   
 84 -  def getUpperBound2D(self, i, j): 
 85      if (i > j): 
 86        j,i = i,j 
 87      return self._boundsMat2D[i, j] 
 88   
 89 -  def getLowerBound2D(self, i, j): 
 90      if (j > i): 
 91        j,i = i,j 
 92      return self._boundsMat2D[i,j] 
 93   
 94   
 95 -  def setUpperBound2D(self, i, j, val, checkBounds=False): 
 96      if (checkBounds): self._checkBounds(i,j) 
 97      if (i > j): 
 98        j,i = i,j 
 99      self._boundsMat2D[i,j] = val 
100   
101 -  def setLowerBound2D(self, i, j, val, checkBounds=False): 
102      if (checkBounds): self._checkBounds(i,j) 
103      if (j > i): 
104        j,i = i,j 
105      self._boundsMat2D[i,j] = val 
106   
107 -  def __str__(self): 
108      res = ' '*13 
109      for i,iFeat in enumerate(self._feats): 
110        res += '% 12s '%iFeat.GetFamily() 
111      res += '\n' 
112      for i,iFeat in enumerate(self._feats): 
113        res += '% 12s '%iFeat.GetFamily() 
114        for j,jFeat in enumerate(self._feats): 
115          if j<i: 
116            res += '% 12.3f '%self.getLowerBound(i,j) 
117          elif j>i: 
118            res += '% 12.3f '%self.getUpperBound(i,j) 
119          else: 
120            res += '% 12.3f '%0.0 
121        res += '\n' 
122      return res 
123         
124 -class ExplicitPharmacophore: 
125    """ this is a pharmacophore with explicit point locations and radii 
126    """ 
127 -  def __init__(self, feats=None,radii=None): 
128      if feats and radii: 
129        self._initializeFeats(feats,radii) 
130   
131 -  def _initializeFeats(self,feats,radii): 
132      if len(feats)!=len(radii): 
133        raise ValueError,'len(feats)!=len(radii)' 
134      self._feats = [] 
135      self._radii = [] 
136      for feat,rad in zip(feats,radii): 
137        if isinstance(feat,ChemicalFeatures.MolChemicalFeature): 
138          pos = feat.GetPos() 
139          newFeat = ChemicalFeatures.FreeChemicalFeature(feat.GetFamily(),feat.GetType(), 
140                                             Geometry.Point3D(pos[0],pos[1],pos[2])) 
141        else: 
142          newFeat = feat 
143        self._feats.append(newFeat) 
144        self._radii.append(rad) 
145         
146 -  def getFeatures(self): 
147      return self._feats 
148   
149 -  def getRadii(self): 
150      return self._radii 
151   
152 -  def getFeature(self, i): 
153      return self._feats[i] 
154   
155 -  def getRadius(self, i): 
156      return self._radii[i] 
157   
158 -  def setRadius(self,i,rad): 
159      self._radii[i]=rad 
160   
161 -  def initFromString(self,text): 
162      lines = text.split(r'\n') 
163      self.initFromLines(lines) 
164   
165 -  def initFromFile(self,inF): 
166      self.initFromLines(inF.readlines()) 
167           
168 -  def initFromLines(self,lines): 
169      from Chem import ChemicalFeatures 
170      logger = logging.logger() 
171   
172      import re 
173      spaces = re.compile('[\ \t]+') 
174   
175      feats = [] 
176      rads = [] 
177      for lineNum,line in enumerate(lines): 
178        txt = line.split('#')[0].strip() 
179        if txt: 
180          splitL = spaces.split(txt) 
181          if len(splitL)<5: 
182            logger.error('Input line %d only contains %d fields, 5 are required. Read failed.'%(lineNum,len(splitL))) 
183            return 
184          fName = splitL[0] 
185          try: 
186            xP = float(splitL[1]) 
187            yP = float(splitL[2]) 
188            zP = float(splitL[3]) 
189            rad = float(splitL[4]) 
190          except ValueError: 
191            logger.error('Error parsing a number of line %d. Read failed.'%(lineNum)) 
192            return 
193          feats.append(ChemicalFeatures.FreeChemicalFeature(fName,fName,Geometry.Point3D(xP,yP,zP))) 
194          rads.append(rad) 
195      self._initializeFeats(feats,rads) 
196   
197 -  def __str__(self): 
198      res = '' 
199      for feat,rad in zip(self._feats,self._radii): 
200        res += '% 12s '%feat.GetFamily() 
201        p = feat.GetPos() 
202        res += '   % 8.4f % 8.4f % 8.4f    '%(p.x,p.y,p.z) 
203        res += '% 5.2f'%rad 
204        res += '\n' 
205      return res 
206