Package rdkit :: Package Chem :: Package Pharm3D :: Module Pharmacophore

[hide private]
[frames] | no frames]

# Source Code for Module rdkit.Chem.Pharm3D.Pharmacophore

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