sidechain_check script

sidechain_check.py (text/x-python), 13.09 KB, created by Jeffrey Gardner (RETIRED) on 2006-03-08 10:21:29 UTC

from pymol.wizard import Wizard
from pymol import cmd,stored
from pymol.parsing import QuietException
import sys

default_hbond_dist = 3.2

class Sidechain_check(Wizard):

   def __init__(self):

      Wizard.__init__(self)
      self.restart()
      cmd.extend('mgz',mgzoom)
      cmd.extend('mgstart',mgstart)
      cmd.extend('flip',flip)
      cmd.extend('chp',cycle_his_protonation)
      cmd.extend('next',next)
      cmd.extend('add_h_bond',_add_naive_hbonds)
      self.surface_state = 0
      

# generic set routines

   def get_panel(self):

      hbond_vals = [self.hbond_dist + i for i in (-3,-2,-1.5,-1.4,-1.3,-1.2,-1.0,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0,
                                                  0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,2,3)]

      self.menu['hbond_dist'] = [[1,str(v),'cmd.get_wizard().set_hbond_dist(%f)'%v] for v in hbond_vals]

      
      if self.surface_state:
         surface_label = 'Hide Surface'
      else:
         surface_label = 'Show Surface'
      return [
         [ 1, 'Check Side-chains',''],
         [ 3, 'H-Bond Cutoff (%s)'%self.hbond_dist,'hbond_dist'],
         [ 2, 'Next Residue' , 'cmd.get_wizard().next()'],
         [ 2, 'Flip', 'cmd.get_wizard().flip()'],
         [ 2, 'Cycle His' , 'cmd.get_wizard().cycle_his_protonation()'],
         [ 2, surface_label , 'cmd.get_wizard().toggle_surface()'],
         [ 2, 'Restart' , 'cmd.get_wizard().restart()'],
         [ 2, 'Done','cmd.set_wizard()'],
         ]

   def set_hbond_dist(self,val):
      self.hbond_dist = val
      res = None
      try:
         res = stored.curr_res
      except AttributeError:
         # no res yet
         pass
      if res is not None:
         add_naive_hbonds(stored.curr_res,hbond_dist=self.hbond_dist)
         
      cmd.refresh_wizard()

   def toggle_surface(self):
      cmd.set('transparency',0.65)
      if self.surface_state:
         cmd.hide('surface')
         self.surface_state = 0
      else:
         cmd.show('surface')
         self.surface_state = 1
      cmd.refresh_wizard()
   
   def restart(self):
      self.cleanup()
      mgstart()
      cmd.show()
      cmd.zoom()
      cmd.refresh_wizard()
   
   def next(self):
      next(self.surface_state,hbond_dist=self.hbond_dist)
      cmd.refresh_wizard()

   def flip(self):
      flip(hbond_dist=self.hbond_dist)
      cmd.refresh_wizard()

   def cycle_his_protonation(self):
      cycle_his_protonation(hbond_dist=self.hbond_dist)
   
   def cleanup(self):
      self.hbond_dist = default_hbond_dist
      self.clear()
      
   def clear(self):
      cmd.delete("HBA")
      cmd.delete("HBD")
      cmd.delete("don")
      cmd.delete("acc")
      cmd.hide("sticks")
      cmd.delete("curr_res")
      self.prompt = None
      try:
         del stored.curr_res
         del stored.tmp
      except AttributeError:
         pass
      
   def get_prompt(self):
      try:
         if type(stored.curr_res) == type(''):
            self.prompt = [ 'Editing residue ' + stored.curr_res ]
         else:
            self.prompt = [ 'Please click "Next Residue" to begin...',
                            'Note that you may need to add Hydrogens',
                            'with "h_add" before you start']
      except AttributeError:
         self.prompt = [ 'Please click "Next Residue" to begin...',
                            'Note that you may need to add Hydrogens',
                            'with "h_add" before you start']
      return self.prompt

###############################################
##                                           ##
## End of class methods                      ##
##                                           ##
###############################################

def mgzoom(i = None, show_surface = None,hbond_dist = default_hbond_dist):
    if i is None:
        try:
            i = stored.curr_res
        except:
            pass
    if type(i) == type(1):
        i = str(i)
    stored.curr_res = i
    print "Current residue: %s" % (i)
    cmd.hide()
    cmd.color('green','(elem C)')
    #
    # Notes:
    #
    # The 'byres' is necessary to show the complete residues.
    # Otherwise we'd just see the nearby atoms.
    #
    # We need to show i and nearby residues as lines before we
    # do anything else.  This makes sure that pymol will actually
    # draw the peptide bonds.  Otherwise, they are sometimes left
    # out.
    #

    visualization_buffer = max(5,hbond_dist+2)
    cmd.show('lines','(' + i + '/)')
    cmd.show('lines','(byres (' + i + '/ around %s))'%visualization_buffer)
    cmd.show('spheres','((r. hoh) and (' + i + '/ around %s))'%visualization_buffer) # water
    cmd.show('sticks','(byres ((resn asn,gln,his,hie,hid,hip) within %s of ('%visualization_buffer + i + '/)))')
    # also show nearby waters
    cmd.show('nb_spheres','(resn hoh within %s of ('%visualization_buffer + i + '/))')
    cmd.zoom('(byres (' + i + '/ around %s))'%(visualization_buffer+3))
    cmd.color('purple','(' + i + '/ and elem C)')
    #
    # Select the correct bond for torsion flipping
    #
    n = get_residue_type(i)
    if n == "GLN":
        cmd.edit('(' + i + '/CD)', '(' + i + '/CG)')
    elif n == "ASN":
        cmd.edit('(' + i + '/CG)', '(' + i + '/CB)')
    elif n  in ('HIS','HID','HIE','HIP'):
        cmd.edit('(' + i + '/CG)', '(' + i + '/CB)')
    add_naive_hbonds(i,hbond_dist=hbond_dist)
    if show_surface:
       cmd.show('surface')

def add_naive_hbonds(i,only_pertinent_bonds=1,hbond_dist=default_hbond_dist):
    n = get_residue_type(i)

    #
    # sel_map maps residue types to selection strings.
    # it tells us the atoms for which we will care about hydrogen bonds.
    # these expect to have two strings substituted in, both the same resi.
    sel_map = {'GLN':'(%s/OE1,NE2 or (hydro and neighbor %s/OE1,NE2))',
               'ASN':'(%s/OD1,ND2 or (hydro and neighbor %s/OD1,ND2))',
               'HIS':'(%s/ND1,NE2 or (hydro and neighbor %s/ND1,NE2))',
               'HIE':'(%s/ND1,NE2 or (hydro and neighbor %s/ND1,NE2))',
               'HID':'(%s/ND1,NE2 or (hydro and neighbor %s/ND1,NE2))',
               'HIP':'(%s/ND1,NE2 or (hydro and neighbor %s/ND1,NE2))',
               }

    if only_pertinent_bonds and n in sel_map.keys():
        return _add_naive_hbonds(i,sel_map[n] % (i,i),hbond_dist=hbond_dist)
    else:
        return _add_naive_hbonds(i,hbond_dist=hbond_dist)


def _add_naive_hbonds(i,sel = '(all)',hbond_dist = default_hbond_dist):
    """\
    Add Hydrogen bonds to i.

    another method of adding naive hydrogen bonds.
    the basic idea for this was stolen from
    http://www.rubor.de/bioinf/pml/hbond2.pml
    which (i assume) came from something warren delano posted on the pymol
    mailing list.
    NOTE: we use 3.2 as the distance cutoff for hydrogen bonding.  The
    distance between the Hydrogen and the heavy atom should really be
    cutoff at 2.6, but we'll use the average heavy atom distance instead,
    since you can't always trust hydrogen positions.
    """
    if type(i) != type('i'):
        i = str(i)
    #cmd.select('don', '(elem n,o and (neighbor hydro)) or (resn hoh) or (resn wat)')
    cmd.select('don', 'elem n,o and (neighbor hydro)')
    #cmd.select('don', 'hydro and neighbor (elem n,o)')
    #cmd.select('acc', 'elem o or (elem n and not (neighbor hydro)) or (resn hoh) or (resn wat)')
    cmd.select('acc', 'elem o or (elem n and not (neighbor hydro))')
    cmd.select('curr_res',sel)
    #cmd.hide('spheres','not resn hoh')
    #cmd.show('spheres','don or acc within %f of curr_res'%hbond_dist)
    try:
        cmd.dist('HBA',
                 '%s/ and acc and %s' % (i,sel),
                 'don and not %s/' % (i),
                 hbond_dist)
        cmd.show('dashes','HBA')
        cmd.show('labels','HBA')
    except QuietException:
        #
        # This happens when there's nothing in one of the selections.
        #
        pass
    try:
        cmd.dist('HBD',
                 '%s/ and don and %s' % (i,sel),
                 'acc and not %s/' % (i),
                 hbond_dist)
        cmd.show('dashes','HBD')
        cmd.show('labels','HBD')
    except QuietException:
        #
        # This happend when there's nothing in one of the selections.
        #
        pass

    #cmd.hide('(hydro)')
    #cmd.delete('don')
    #cmd.delete('acc')

def get_residue_type(i):
    #
    # There *MUST* be a better way ...
    #
    if type(i) == type(1):
        i = str(i)
    stored.tmp = None
    cmd.iterate('(' + i + '/)','stored.tmp = resn')
    return stored.tmp

def next(show_surface = None,hbond_dist=default_hbond_dist):
   res = stored.interesting_residues.pop()
   mgzoom(res, show_surface,
          hbond_dist=hbond_dist)
   #
   # cheap hack to make sure that we're displaying the
   # right version of HIS, HID, HIE or HIP ..
   # if we're listed as HIS, we'll force it to be HIP.
   # that requires cycling through the protonation
   # 3 times HIS->HIE->HID->HIP
   # this happens to be necessary if you use h_add to
   # set up histidine protonation.
   #
   if get_residue_type(res) in ('HIS',):
      cycle_his_protonation()
      cycle_his_protonation()
      cycle_his_protonation()
      pass
    

def flip(delete_old_hbonds = 1,hbond_dist=default_hbond_dist):
   cmd.delete('HBA')
   cmd.delete('HBD')
   res_type = get_residue_type(stored.curr_res)
   if res_type in ('ASN','GLN','HIS','HID','HIE','HIP'):
      cmd.torsion(180)
   add_naive_hbonds(stored.curr_res,hbond_dist=hbond_dist)

def mgstart(addhydrogen = None):
    """
    You should use protonate to add hydrogens.  Failing that, set
    addhydrogen to something non-false and PyMOL will try to do it for
    you.
    """
    cmd.set("sphere_scale",0.25)
    if addhydrogen:
        cmd.h_add()
    stored.list = []
    cmd.iterate("(resn asn,gln,his,hie,hid,hip)","stored.list.append((resi,resn))")
    stored.interesting_residues = {}
    for item in stored.list:
        i,n = item
        stored.interesting_residues[i] = 1
    stored.interesting_residues = stored.interesting_residues.keys()
    stored.interesting_residues.sort()
    stored.interesting_residues.reverse()


def cycle_his_protonation(i = None,hbond_dist=default_hbond_dist):
    """Flip between HIE, HID and HIP

    change the protonation state of a his residue.
    HIE : NE2 protonated
    HID : ND1 protonated
    HIP : both protonated

    We'll cycle from HIE --> HID --> HIP on successive calls.
    """
    if i is None:
        try:
            i = stored.curr_res
        except:
            pass
    if type(i) == type(1):
        i = str(i)
    curr_state = get_residue_type(i)
    #
    # There are a few things we need to do here:
    #
    # 1) Change the name from HIE to HID
    # 2) Change the valences
    # 3) Add the Hydrogens
    #
    # Location of double bonds:
    #  - All three have a 2x bond between CG and CD2
    #  - HIE has a 2x bond between ND1 and CE1
    #  - HID has a 2x bond between NE2 and CE1
    #
    # Note: PyMOL has no set_valence command, so I have to use cycle_valence
    # and trust that nobody has screwed with the valences.
    # This is sub-optimal.
    #
    # Note: Do not forget to h_fill on the bonds that 'stay the same' .. they
    # stay as single bonds, but the protonation state changes.
    #
    # Note: for HIP, we put a +1 charge on the ND1 Nitrogen.  Don't forget to
    # remove it when you cycle to HIE.
    #
    cmd.remove('(hydro and neighbor %s/ND1,NE2)' % i)
    if curr_state == 'HIE':
        print "HIE --> HID"
        cmd.alter('%s/' % i,'resn="HID"')
        # ND1=CE1 must change from 2x to 1x.  2x->3x->1x == two cycles
        cmd.edit('(' + i + '/ND1)', '(' + i + '/CE1)')
        cmd.cycle_valence()
        cmd.cycle_valence()
        # NE2-CE1 must change from 1x to 2x.  1x->2x == one cycle
        cmd.edit('(' + i + '/NE2)', '(' + i + '/CE1)')
        cmd.cycle_valence()
    elif curr_state == 'HID':
        print "HID --> HIP"
        cmd.alter('%s/' % i,'resn="HIP"')
        # NE2-CE1 must change from 2x to 1x.  2x->3x->1x == two cycles
        cmd.edit('(' + i + '/NE2)', '(' + i + '/CE1)')
        cmd.cycle_valence()
        cmd.cycle_valence()
        # ND1-CE1 must change from 1x to 2x.  1x->2x == one cycle
        cmd.edit('(' + i + '/ND1)', '(' + i + '/CE1)')
        cmd.cycle_valence()
        #
        # Add the +1 charge to ND1 as mentioned above
        #
        cmd.alter('(' + i + '/ND1)','formal_charge=1')
    elif curr_state == 'HIP' or curr_state == 'HIS':
        if curr_state == 'HIS':
           print "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
           print "!!!!! Assuming HIS means HIP !!!!!!"
           print "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
        print "HIP --> HIE"
        cmd.alter('%s/' % i,'resn="HIE"')
        #
        # Remove the +1 charge from ND1 as mentioned above
        #
        cmd.alter('(' + i + '/ND1)','formal_charge=0')

    # Make sure we've h_fill'd the appropriate bonds.
    cmd.sort()
    cmd.edit('(' + i + '/ND1)', '(' + i + '/CE1)')
    cmd.h_fill()
    cmd.edit('(' + i + '/NE2)', '(' + i + '/CE1)')
    cmd.h_fill()
    add_naive_hbonds(stored.curr_res,hbond_dist=hbond_dist)
    #
    # When we're all done, we want the CG-CB bond selected so that 'flip'
    # will still work.
    #
    cmd.edit('(' + i + '/CG)', '(' + i + '/CB)')