Source Code for Module rdkit.Chem.FeatFinderCLI

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