Source Code for Module rdkit.Chem.Pharm3D.Pharmacophore

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