Source Code for Module rdkit.Chem.Pharm2D.Matcher

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