Source Code for Module rdkit.Chem.FeatFinderCLI

  1  # $Id: FeatFinderCLI.py 997 2009-02-25 06:12:43Z glandrum $ 
  2  # 
  3  #  Copyright (C) 2005-2006 Rational Discovery LLC 
  4  # 
  5  #   @@ All Rights Reserved  @@ 
  6  # 
  7  _version = "$Rev: 997 $" 
  8  _splashMessage=""" 
  9  -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* 
 10    FeatFinderCLI version %s 
 11   
 12    Copyright (C) 2005 Rational Discovery LLC 
 13   
 14    This software is copyrighted.  The software may not be copied, 
 15    reproduced, translated or reduced to any electronic medium or 
 16    machine-readable form without the prior written consent of 
 17    Rational Discovery LLC. 
 18  -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* 
 19  """%_version 
 20  from rdkit import Chem 
 21  from rdkit.Chem import ChemicalFeatures 
 22  from rdkit import RDLogger 
 23  logger = RDLogger.logger() 
 24  import sys,os 
 25  import re 
 26  splitExpr = re.compile(r'[ \t,]') 
 27   
 28 -def GetAtomFeatInfo(factory,mol): 
 29    res = [None]*mol.GetNumAtoms() 
 30    feats = factory.GetFeaturesForMol(mol) 
 31    for feat in feats: 
 32      ids = feat.GetAtomIds() 
 33      for id in ids: 
 34        if res[id] is None: 
 35          res[id] = [] 
 36        res[id].append("%s-%s"%(feat.GetFamily(),feat.GetType())) 
 37    return res 
 38   
 39  if __name__ == '__main__': 
 40 -  def Usage(): 
 41      message=""" 
 42      Usage: FeatFinderCLI [-r] <fdefFilename> <smilesFilename> 
 43       
 44      NOTE: 
 45        - the smiles file should have SMILES in the first column 
 46   
 47      """ 
 48      print >>sys.stderr,message 
 49   
 50   
 51    import getopt 
 52    args,extras = getopt.getopt(sys.argv[1:],'r') 
 53    reverseIt=False 
 54    for arg,val in args: 
 55      if arg=='-r': 
 56        reverseIt=True 
 57         
 58    if len(extras)<2: 
 59      Usage() 
 60      sys.exit(-1) 
 61    print >>sys.stderr,_splashMessage 
 62    fdefFilename = extras[0] 
 63    if not os.path.exists(fdefFilename): 
 64      logger.error("Fdef file %s does not exist."%fdefFilename) 
 65      sys.exit(-1) 
 66    try: 
 67      factory = ChemicalFeatures.BuildFeatureFactory(fdefFilename) 
 68    except: 
 69      logger.error("Could not parse Fdef file %s."%fdefFilename,exc_info=True) 
 70      sys.exit(-1) 
 71       
 72    smilesFilename = extras[1] 
 73    if not os.path.exists(smilesFilename): 
 74      logger.error("Smiles file %s does not exist."%smilesFilename) 
 75      sys.exit(-1) 
 76   
 77    try: 
 78      inF = file(smilesFilename,'r') 
 79    except: 
 80      logger.error("Could not open smiles file %s."%smilesFilename,exc_info=True) 
 81      sys.exit(-1) 
 82   
 83    lineNo=0 
 84    for line in inF.readlines(): 
 85      lineNo+=1 
 86      line = line.strip() 
 87      smi = splitExpr.split(line)[0].strip() 
 88      try: 
 89        mol = Chem.MolFromSmiles(smi) 
 90      except: 
 91        mol = None 
 92   
 93      if mol: 
 94        print 'Mol-%d\t%s'%(lineNo,smi) 
 95   
 96        if not reverseIt: 
 97          featInfo = GetAtomFeatInfo(factory,mol) 
 98          for i,v in enumerate(featInfo): 
 99            print '\t% 2s(%d)'%(mol.GetAtomWithIdx(i).GetSymbol(),i+1), 
100            if v: 
101              print '\t',', '.join(v) 
102            else: 
103              print 
104        else: 
105          feats = factory.GetFeaturesForMol(mol) 
106          for feat in feats: 
107            print '\t%s-%s: '%(feat.GetFamily(),feat.GetType()), 
108            print ', '.join([str(x) for x in feat.GetAtomIds()]) 
109      else: 
110        logger.warning("Could not process smiles '%s' on line %d."%(smi,lineNo)) 
111