Package rdkit :: Package Chem :: Package Pharm2D :: Module Matcher

[hide private]
[frames] | no frames]

# Source Code for Module rdkit.Chem.Pharm2D.Matcher

      1  # $Id$
      2  #
      3  # Copyright (C) 2003-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  """ functionality for finding pharmacophore matches in molecules
     12
     13
     14    See Docs/Chem/Pharm2D.triangles.jpg for an illustration of the way
     15    pharmacophores are broken into triangles and labelled.
     16
     17    See Docs/Chem/Pharm2D.signatures.jpg for an illustration of bit
     18    numbering
     19
     20  """
     21  from rdkit import Chem
     22  from rdkit.Chem.Pharm2D import Utils
     23
     24  import types
     25
     26 -class MatchError(Exception):
     27    pass
     28
     29  _verbose = 0
     30 -def GetAtomsMatchingBit(sigFactory,bitIdx,mol,dMat=None,justOne=0,matchingAtoms=None):
     31    """ Returns a list of lists of atom indices for a bit
     32
     33      **Arguments**
     34
     35        - sigFactory: a SigFactory
     36
     37        - bitIdx: the bit to be queried
     38
     39        - mol: the molecule to be examined
     40
     41        - dMat: (optional) the distance matrix of the molecule
     42
     43        - justOne: (optional) if this is nonzero, only the first match
     44          will be returned.
     45
     46        - matchingAtoms: (optional) if this is nonzero, it should
     47          contain a sequence of sequences with the indices of atoms in
     48          the molecule which match each of the patterns used by the
     49          signature.
     50
     51      **Returns**
     52
     53        a list of tuples with the matching atoms
     54    """
     55    assert sigFactory.shortestPathsOnly,'not implemented for non-shortest path signatures'
     56    nPts,featCombo,scaffold = sigFactory.GetBitInfo(bitIdx)
     57    if _verbose:
     58      print('info:',nPts)
     59      print('\t',featCombo)
     60      print('\t',scaffold)
     61
     62    if matchingAtoms is None:
     63      matchingAtoms = sigFactory.GetMolFeats(mol)
     64
     65    # find the atoms that match each features
     66    fams = sigFactory.GetFeatFamilies()
     67    choices = []
     68    for featIdx in featCombo:
     69      tmp = matchingAtoms[featIdx]
     70      if tmp:
     71        choices.append(tmp)
     72      else:
     73        # one of the patterns didn't find a match, we
     74        #  can return now
     75        if _verbose: print('no match found for feature:',featIdx)
     76        return []
     77
     78    if _verbose:
     79      print('choices:')
     80      print(choices)
     81
     82    if dMat is None:
     83      dMat = Chem.GetDistanceMatrix(mol,sigFactory.includeBondOrder)
     84
     85    matches = []
     86    distsToCheck = Utils.nPointDistDict[nPts]
     87
     88    protoPharmacophores = Utils.GetAllCombinations(choices,noDups=1)
     89
     90    res = []
     91    for protoPharm in protoPharmacophores:
     92      if _verbose: print('protoPharm:',protoPharm)
     93      for i in range(len(distsToCheck)):
     94        dLow,dHigh = sigFactory.GetBins()[scaffold[i]]
     95        a1,a2 = distsToCheck[i]
     96        #
     97        # FIX: this is making all kinds of assumptions about
     98        #  things being single-atom matches (or at least that
     99        #  only the first atom matters
    100        #
    101        idx1,idx2 = protoPharm[a1][0],protoPharm[a2][0]
    102        dist = dMat[idx1,idx2]
    103        if _verbose: print('\t dist: %d->%d = %d (%d,%d)'%(idx1,idx2,dist,dLow,dHigh))
    104        if dist < dLow or dist >= dHigh:
    105          break
    106      else:
    107        if _verbose: print('Found one')
    108        # we found it
    109        protoPharm.sort()
    110        protoPharm = tuple(protoPharm)
    111        if protoPharm not in res:
    112          res.append(protoPharm)
    113          if justOne: break
    114    return res
    115
    116  if __name__ == '__main__':
    117    from rdkit import Chem
    118    from rdkit.Chem.Pharm2D import SigFactory,Generate
    119
    120    factory = SigFactory.SigFactory()
    121    factory.SetBins([(1,2),(2,5),(5,8)])
    122    factory.SetPatternsFromSmarts(['O','N'])
    123    factory.SetMinCount(2)
    124    factory.SetMaxCount(3)
    125    sig = factory.GetSignature()
    126
    127    mol = Chem.MolFromSmiles('OCC(=O)CCCN')
    128    Generate.Gen2DFingerprint(mol,sig)
    129    print('onbits:',list(sig.GetOnBits()))
    130
    131    _verbose=0
    132    for bit in sig.GetOnBits():
    133      atoms = GetAtomsMatchingBit(sig,bit,mol)
    134      print('\tBit %d: '%(bit),atoms)
    135
    136
    137    print('--------------------------')
    138    sig = factory.GetSignature()
    139    sig.SetIncludeBondOrder(1)
    140    Generate.Gen2DFingerprint(mol,sig)
    141    print('onbits:',list(sig.GetOnBits()))
    142
    143    for bit in sig.GetOnBits():
    144      atoms = GetAtomsMatchingBit(sig,bit,mol)
    145      print('\tBit %d: '%(bit),atoms)
    146