Package rdkit :: Package Chem :: Package Features :: Module ShowFeats

[hide private]
[frames] | no frames]

# Source Code for Module rdkit.Chem.Features.ShowFeats

      1  # $Id: ShowFeats.py 537 2007-08-20 14:54:35Z landrgr1 $
      2  #
      3  # Created by Greg Landrum Aug 2006
      4  #
      5  #
      6  from __future__ import print_function
      7
      8  _version = "0.3.2"
      9
     10
     11  _usage="""
     12     ShowFeats [optional args] <filenames>
     13
     14    if "-" is provided as a filename, data will be read from stdin (the console)
     15  """
     16
     17  _welcomeMessage="This is ShowFeats version %s"%(_version)
     18
     19
     20  import math
     21  #set up the logger:
     22  from rdkit import RDLogger as logging
     23  logger = logging.logger()
     24  logger.setLevel(logging.INFO)
     25
     26  from rdkit import Geometry
     27  from rdkit.Chem.Features import FeatDirUtilsRD as FeatDirUtils
     28
     29  _featColors = {
     30    'Donor':(0,1,1),
     31    'Acceptor':(1,0,1),
     32    'NegIonizable':(1,0,0),
     33    'PosIonizable':(0,0,1),
     34    'ZnBinder':(1,.5,.5),
     35    'Aromatic':(1,.8,.2),
     36    'LumpedHydrophobe':(.5,.25,0),
     37    'Hydrophobe':(.5,.25,0),
     38    }
     39
     40 -def _getVectNormal(v,tol=1e-4):
     41    if math.fabs(v.x)>tol:
     42      res = Geometry.Point3D(v.y,-v.x,0)
     43    elif math.fabs(v.y)>tol:
     44      res = Geometry.Point3D(-v.y,v.x,0)
     45    elif math.fabs(v.z)>tol:
     46      res = Geometry.Point3D(1,0,0)
     47    else:
     48      raise ValueError('cannot find normal to the null vector')
     49    res.Normalize()
     50    return res
     51
     52  _canonArrowhead=None
     53 -def _buildCanonArrowhead(headFrac,nSteps,aspect):
     54    global _canonArrowhead
     55    startP = RDGeometry.Point3D(0,0,headFrac)
     56    _canonArrowhead=[startP]
     57
     58    scale = headFrac*aspect
     59    baseV = RDGeometry.Point3D(scale,0,0)
     60    _canonArrowhead.append(baseV)
     61
     62    twopi = 2*math.pi
     63    for i in range(1,nSteps):
     64      v = RDGeometry.Point3D(scale*math.cos(i*twopi),scale*math.sin(i*twopi),0)
     65      _canonArrowhead.append(v)
     66
     67
     68  _globalArrowCGO=[]
     69  _globalSphereCGO=[]
     70  # taken from pymol's cgo.py
     71  BEGIN=2
     72  END=3
     73  TRIANGLE_FAN=6
     74  COLOR=6
     75  VERTEX=4
     76  NORMAL=5
     77  SPHERE=7
     78  CYLINDER=9
     79  ALPHA=25
     80
     81 -def _cgoArrowhead(viewer,tail,head,radius,color,label,headFrac=0.3,nSteps=10,aspect=.5):
     82    global _globalArrowCGO
     83    delta = head-tail
     84    normal = _getVectNormal(delta)
     85    delta.Normalize()
     86
     87    dv = head-tail
     88    dv.Normalize()
     89    dv *= headFrac
     90    startP = head
     91
     92    normal*=headFrac*aspect
     93
     94    cgo = [BEGIN,TRIANGLE_FAN,
     95           COLOR,color[0],color[1],color[2],
     96            NORMAL,dv.x,dv.y,dv.z,
     97           VERTEX,head.x+dv.x,head.y+dv.y,head.z+dv.z]
     98    base = [BEGIN,TRIANGLE_FAN,
     99            COLOR,color[0],color[1],color[2],
    100            NORMAL,-dv.x,-dv.y,-dv.z,
    101            VERTEX,head.x,head.y,head.z]
    102    v = startP+normal
    103    cgo.extend([NORMAL,normal.x,normal.y,normal.z])
    104    cgo.extend([VERTEX,v.x,v.y,v.z])
    105    base.extend([VERTEX,v.x,v.y,v.z])
    106    for i in range(1,nSteps):
    107      v = FeatDirUtils.ArbAxisRotation(360./nSteps*i,delta,normal)
    108      cgo.extend([NORMAL,v.x,v.y,v.z])
    109      v += startP
    110      cgo.extend([VERTEX,v.x,v.y,v.z])
    111      base.extend([VERTEX,v.x,v.y,v.z])
    112
    113    cgo.extend([NORMAL,normal.x,normal.y,normal.z])
    114    cgo.extend([VERTEX,startP.x+normal.x,startP.y+normal.y,startP.z+normal.z])
    115    base.extend([VERTEX,startP.x+normal.x,startP.y+normal.y,startP.z+normal.z])
    116    cgo.append(END)
    117    base.append(END)
    118    cgo.extend(base)
    119
    120    #viewer.server.renderCGO(cgo,label)
    121    _globalArrowCGO.extend(cgo)
    122
    123 -def ShowArrow(viewer,tail,head,radius,color,label,transparency=0,includeArrowhead=True):
    124    global _globalArrowCGO
    125    if transparency:
    126      _globalArrowCGO.extend([ALPHA,1-transparency])
    127    else:
    128      _globalArrowCGO.extend([ALPHA,1])
    129    _globalArrowCGO.extend([CYLINDER,tail.x,tail.y,tail.z,
    130                            head.x,head.y,head.z,
    131                            radius*.10,
    132                            color[0],color[1],color[2],
    133                            color[0],color[1],color[2],
    134                            ])
    135
    136    if includeArrowhead:
    137      _cgoArrowhead(viewer,tail,head,radius,color,label)
    138
    139
    140 -def ShowMolFeats(mol,factory,viewer,radius=0.5,confId=-1,showOnly=True,
    141                   name='',transparency=0.0,colors=None,excludeTypes=[],
    142                   useFeatDirs=True,featLabel=None,dirLabel=None,includeArrowheads=True,
    143                   writeFeats=False,showMol=True,featMapFile=False):
    144    global _globalSphereCGO
    145    if not name:
    146      if mol.HasProp('_Name'):
    147        name =  mol.GetProp('_Name')
    148      else:
    149        name = 'molecule'
    150    if not colors:
    151      colors = _featColors
    152
    153    if showMol:
    154      viewer.ShowMol(mol,name=name,showOnly=showOnly,confId=confId)
    155
    156    molFeats=factory.GetFeaturesForMol(mol)
    157    if not featLabel:
    158      featLabel='%s-feats'%name
    159      viewer.server.resetCGO(featLabel)
    160    if not dirLabel:
    161      dirLabel=featLabel+"-dirs"
    162      viewer.server.resetCGO(dirLabel)
    163
    164    for i,feat in enumerate(molFeats):
    165      family=feat.GetFamily()
    166      if family in excludeTypes:
    167        continue
    168      pos = feat.GetPos(confId)
    169      color = colors.get(family,(.5,.5,.5))
    170      nm = '%s(%d)'%(family,i+1)
    171
    172      if transparency:
    173        _globalSphereCGO.extend([ALPHA,1-transparency])
    174      else:
    175        _globalSphereCGO.extend([ALPHA,1])
    176      _globalSphereCGO.extend([COLOR,color[0],color[1],color[2],
    177                               SPHERE,pos.x,pos.y,pos.z,
    178                               radius])
    179      if writeFeats:
    180        aidText = ' '.join([str(x+1) for x in feat.GetAtomIds()])
    181        print('%s\t%.3f\t%.3f\t%.3f\t1.0\t# %s'%(family,pos.x,pos.y,pos.z,aidText))
    182
    183      if featMapFile:
    184        print("  family=%s pos=(%.3f,%.3f,%.3f) weight=1.0"%(family,pos.x,pos.y,pos.z),end='',file=featMapFile)
    185
    186      if useFeatDirs:
    187        ps = []
    188        if family=='Aromatic':
    189          ps,fType = FeatDirUtils.GetAromaticFeatVects(mol.GetConformer(confId),
    190                                                                 feat.GetAtomIds(),pos,
    191                                                                 scale=1.0)
    192        elif family=='Donor':
    193          aids = feat.GetAtomIds()
    194          if len(aids)==1:
    195            featAtom=mol.GetAtomWithIdx(aids[0])
    196            hvyNbrs=[x for x in featAtom.GetNeighbors() if x.GetAtomicNum()!=1]
    197            if len(hvyNbrs)==1:
    198              ps,fType = FeatDirUtils.GetDonor1FeatVects(mol.GetConformer(confId),
    199                                                         aids,scale=1.0)
    200            elif len(hvyNbrs)==2:
    201              ps,fType = FeatDirUtils.GetDonor2FeatVects(mol.GetConformer(confId),
    202                                                         aids,scale=1.0)
    203            elif len(hvyNbrs)==3:
    204              ps,fType = FeatDirUtils.GetDonor3FeatVects(mol.GetConformer(confId),
    205                                                         aids,scale=1.0)
    206        elif family=='Acceptor':
    207          aids = feat.GetAtomIds()
    208          if len(aids)==1:
    209            featAtom=mol.GetAtomWithIdx(aids[0])
    210            hvyNbrs=[x for x in featAtom.GetNeighbors() if x.GetAtomicNum()!=1]
    211            if len(hvyNbrs)==1:
    212              ps,fType = FeatDirUtils.GetAcceptor1FeatVects(mol.GetConformer(confId),
    213                                                            aids,scale=1.0)
    214            elif len(hvyNbrs)==2:
    215              ps,fType = FeatDirUtils.GetAcceptor2FeatVects(mol.GetConformer(confId),
    216                                                         aids,scale=1.0)
    217            elif len(hvyNbrs)==3:
    218              ps,fType = FeatDirUtils.GetAcceptor3FeatVects(mol.GetConformer(confId),
    219                                                         aids,scale=1.0)
    220
    221        for tail,head in ps:
    222          ShowArrow(viewer,tail,head,radius,color,dirLabel,
    223                    transparency=transparency,includeArrowhead=includeArrowheads)
    224          if featMapFile:
    225            vect = head-tail
    226            print('dir=(%.3f,%.3f,%.3f)'%(vect.x,vect.y,vect.z),end='',file=featMapFile)
    227
    228      if featMapFile:
    229        aidText = ' '.join([str(x+1) for x in feat.GetAtomIds()])
    230        print('# %s'%(aidText),file=featMapFile)
    231
    232
    233
    234
    235  # --- ----  --- ----  --- ----  --- ----  --- ----  --- ----
    236  import sys,os,getopt
    237  from rdkit import RDConfig
    238  from optparse import OptionParser
    239  parser=OptionParser(_usage,version='%prog '+_version)
    240
    241  parser.add_option('-x','--exclude',default='',
    242                    help='provide a list of feature names that should be excluded')
    243  parser.add_option('-f','--fdef',default=os.path.join(RDConfig.RDDataDir,'BaseFeatures.fdef'),
    244                    help='provide the name of the feature definition (fdef) file.')
    245  parser.add_option('--noDirs','--nodirs',dest='useDirs',default=True,action='store_false',
    246                    help='do not draw feature direction indicators')
    247  parser.add_option('--noHeads',dest='includeArrowheads',default=True,action='store_false',
    248                    help='do not draw arrowheads on the feature direction indicators')
    249  parser.add_option('--noClear','--noClear',dest='clearAll',default=False,action='store_true',
    250                    help='do not clear PyMol on startup')
    251  parser.add_option('--noMols','--nomols',default=False,action='store_true',
    252                    help='do not draw the molecules')
    253  parser.add_option('--writeFeats','--write',default=False,action='store_true',
    254                    help='print the feature information to the console')
    255  parser.add_option('--featMapFile','--mapFile',default='',
    256                    help='save a feature map definition to the specified file')
    257  parser.add_option('--verbose',default=False,action='store_true',
    258                    help='be verbose')
    259
    260  if __name__=='__main__':
    261    from rdkit import Chem
    262    from rdkit.Chem import AllChem
    263    from rdkit.Chem.PyMol import MolViewer
    264
    265    options,args = parser.parse_args()
    266    if len(args)<1:
    267      parser.error('please provide either at least one sd or mol file')
    268
    269    try:
    270      v = MolViewer()
    271    except:
    272      logger.error('Unable to connect to PyMol server.\nPlease run ~landrgr1/extern/PyMol/launch.sh to start it.')
    273      sys.exit(1)
    274    if options.clearAll:
    275      v.DeleteAll()
    276
    277
    278    try:
    279      fdef = open(options.fdef,'r').read()
    280    except IOError:
    281      logger.error('ERROR: Could not open fdef file %s'%options.fdef)
    282      sys.exit(1)
    283
    284    factory = AllChem.BuildFeatureFactoryFromString(fdef)
    285
    286    if options.writeFeats:
    287       print('# Family   \tX    \tY   \tZ   \tRadius\t # Atom_ids')
    288
    289    if options.featMapFile:
    290      if options.featMapFile=='-':
    291        options.featMapFile=sys.stdout
    292      else:
    293        options.featMapFile=file(options.featMapFile,'w+')
    294      print('# Feature map generated by ShowFeats v%s'%_version, file=options.featMapFile)
    295      print("ScoreMode=All", file=options.featMapFile)
    296      print("DirScoreMode=Ignore", file=options.featMapFile)
    297      print("BeginParams", file=options.featMapFile)
    298      for family in factory.GetFeatureFamilies():
    299        print("   family=%s width=1.0 radius=3.0"%family, file=options.featMapFile)
    300      print("EndParams", file=options.featMapFile)
    301      print("BeginPoints", file=options.featMapFile)
    302
    303    i = 1
    304    for midx,molN in enumerate(args):
    305      if molN!='-':
    306        featLabel='%s_Feats'%molN
    307      else:
    308        featLabel='Mol%d_Feats'%(midx+1)
    309
    310      v.server.resetCGO(featLabel)
    311      # this is a big of kludgery to work around what seems to be a pymol cgo bug:
    312      v.server.sphere((0,0,0),.01,(1,0,1),featLabel)
    313      dirLabel=featLabel+"-dirs"
    314
    315      v.server.resetCGO(dirLabel)
    316      # this is a big of kludgery to work around what seems to be a pymol cgo bug:
    317      v.server.cylinder((0,0,0),(.01,.01,.01),.01,(1,0,1),dirLabel)
    318
    319      if molN != '-':
    320        try:
    321          ms = Chem.SDMolSupplier(molN)
    322        except:
    323          logger.error('Problems reading input file: %s'%molN)
    324          ms = []
    325      else:
    326        ms = Chem.SDMolSupplier()
    327        ms.SetData(sys.stdin.read())
    328
    329      for m in ms:
    330        nm = 'Mol_%d'%(i)
    331        if m.HasProp('_Name'):
    332          nm += '_'+m.GetProp('_Name')
    333        if options.verbose:
    334          if m.HasProp('_Name'):
    335            print("#Molecule: %s"%m.GetProp('_Name'))
    336          else:
    337            print("#Molecule: %s"%nm)
    338        ShowMolFeats(m,factory,v,transparency=0.25,excludeTypes=options.exclude,name=nm,
    339                     showOnly=False,
    340                     useFeatDirs=options.useDirs,
    341                     featLabel=featLabel,dirLabel=dirLabel,
    342                     includeArrowheads=options.includeArrowheads,
    343                     writeFeats=options.writeFeats,showMol=not options.noMols,
    344                     featMapFile=options.featMapFile)
    345        i += 1
    346        if not i%100:
    347          logger.info("Done %d poses"%i)
    348      if ms:
    349        v.server.renderCGO(_globalSphereCGO,featLabel,1)
    350        if options.useDirs:
    351          v.server.renderCGO(_globalArrowCGO,dirLabel,1)
    352
    353    if options.featMapFile:
    354      print("EndPoints",file=options.featMapFile)
    355
    356
    357    sys.exit(0)
    358