Cutting a "binding site apple"

Two weeks ago I had an idea that just couldn't let me go. It all started when was looking at those beautiful binding pocket images of the glucocorticoid receptor back on Felix's blog. Often when I work with Pymol I ask myself how to show the bounded molecule in the binding site in the most interesting way. Usually it's a lot of work involving changing representation modes and hiding residues that are in the way. The last thing is very time-consuming since there is no command for getting this job done easily. You have to either select residues manually and then hide or delete them. Of course you could just move the two clipping planes using the clip command, but some times it's hard to hide residues that are in front without hiding parts of the bounded molecule.

Ok, now we know the problem I'll show you what solution I've worked out for this. I call it the "cut the binding site apple script solution" (yes, I know it's a lame name, but it really hits the nail on the head). Think of it that way - the binding site is an apple where the ligand is the seed and the surrounding residues is the apple flesh. If we cut the apple in the middle we can easily see the bounded ligand. As you probably guessed on your own the cutting is just a simple complement operation on a set of selected binding site atoms. The tricky part was to get the right cutting plane. When I realized that it's perpendicular to the point of current view vector it was just a matter of writing the code. And here it is - my binding pocket surface generator plugin for Pymol.

bp_surface.py(Toggle Plain Text)

# -*- coding: utf-8 -*-

# -------------------------------------------------------------------------------
# bp_surface.py - Binding pocket surface generator for Pymol version 1.0
# -------------------------------------------------------------------------------

# Copyright (C) 2008 by Lightnir - lightnir@gmail.com

# This script is free software; you can redistribute it and#or modify
# it under the terms of the GNU Library General Public License as
# published by the Free Software Foundation; either version 2 of the
# License, or (at your option) any later version.

# This script is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.

# You should have received a copy of the GNU Library General Public
# License along with this program; if not, write to the
# Free Software Foundation, Inc.,
# 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

from pymol import cmd
from Tkinter import *
import Pmw

BPSFrame = None
selections = []
cv = IntVar()

def __init__(self):
    self.menuBar.addcascademenu('Plugin', 'LightnirsPlugins', 'Lightnirs Plugins',
                                    label='Lightnirs Plugins')
    self.menuBar.addmenuitem('LightnirsPlugins', 'command',
                                 'BPsurface',
                                 label='BP surface',
                                 command = lambda s=self: BPSurface(s))
class BPSurface:
    def __init__(self,app):
        global BPSFrame
        global selections
        global cv
        if BPSFrame==None:
            BPSFrame = Toplevel(width=400, height=300)              # create a root window
            BPSFrame.title("BP surface editor")                     # Set window title
            BPSFrame.resizable(0, 0)                                # Disable window resize
            BPSFrame.geometry('-40+40')                             # Set window placement
            self.get_selections()
            # Create Ligand selection Combobox
            BPSFrame.cb1 = Pmw.ComboBox(BPSFrame,
                label_text = 'Ligand selection:',
                labelpos = 'nw',
                scrolledlist_items = selections
            )
            BPSFrame.cb1.grid(row=0, column=0,columnspan=2)

            # Create binding pocket Combobox
            BPSFrame.cb2 = Pmw.ComboBox(BPSFrame,
                label_text = 'Binding pocket selection:',
                labelpos = 'nw',
                scrolledlist_items = selections
            )
            BPSFrame.cb2.grid(row=1, column=0,columnspan=2)

            # Create Slider
            BPSFrame.slider = Scale(BPSFrame, from_=-20, to=20, tickinterval=10, resolution=0.5, label='Z axis cutoff', orient=HORIZONTAL)
            BPSFrame.slider.grid(row=2, column=0, columnspan=2, sticky=N+W+S+E,)

            # Add Checkbox
            cbx = Checkbutton(BPSFrame, text="Create new object", variable=cv, onvalue="1", offvalue="0")
            cbx.select()
            cbx.grid(row=3,column=0, columnspan=2, sticky=W)

            # Add some buttons
            btn1=Button(BPSFrame, text='Refresh', padx=0, pady=0, command = self.refresh_list)
            btn1.grid(row=4, column=0, sticky=N+W+S+E, padx=0, pady=0,ipadx=0,ipady=0)
            btn2=Button(BPSFrame, text='Clear', padx=0, pady=0, command = lambda r=1,d=1: self.clear_flags(r,d))
            btn2.grid(row=4, column=1, sticky=N+W+S+E, padx=0, pady=0,ipadx=0,ipady=0)
            btn3=Button(BPSFrame, text='Apply', padx=0, pady=0, command = self.ok)
            btn3.grid(row=5, column=0, sticky=N+W+S+E,)
            btn4=Button(BPSFrame, text='Close', padx=0, pady=0, command = self.myHide)
            btn4.grid(row=5, column=1, sticky=N+W+S+E)

            # create callback to prevent window kill
            BPSFrame.protocol("WM_DELETE_WINDOW", self.myHide)
            BPSFrame.mainloop()
        else:
            if BPSFrame.state() == "normal":
                self.myHide()
            elif BPSFrame.state() == "withdrawn":
                self.myShow()

    def myShow(self):
        self.refresh_list()
        BPSFrame.deiconify()

    def myHide(self):
        BPSFrame.withdraw()

    def get_selections(self):
        global selections
        selections = []
        for item in cmd.get_names("selections"):
            if item[0]<>"_":
                selections.append(item)

    def refresh_list(self):
        global selections

        # Save old values
        old = selections
        try:
            old_cb1=BPSFrame.cb1.getvalue()[0]
        except:
            old_cb1=''
        try:
            old_cb2=BPSFrame.cb2.getvalue()[0]
        except:
            old_cb2=''
        self.get_selections()

        # Update selectionlist if needed
        if not(old==selections):
            BPSFrame.cb1.setlist(selections)
            BPSFrame.cb2.setlist(selections)

        # Is ligand selection still available?
        if not(old_cb1 in selections):
            # No - clear the combobox entry
            BPSFrame.cb1.component('entryfield').clear()
        else:
            # Yes - set the old value
            BPSFrame.cb1.selectitem(old_cb1,setentry=1)

        # Is binding pocket selection still available?
        if not(old_cb2 in selections):
            # No - clear the combobox entry
            BPSFrame.cb2.component('entryfield').clear()
        else:
            # Yes - set the old value
            BPSFrame.cb2.selectitem(old_cb2,setentry=1)

    def clear_flags(self,r,d):
        cmd.flag('ignore','all','clear')
        if r==1:
            cmd.rebuild('all')
        if d==1:
            cmd.delete('bp_cutoff')
            cmd.delete('bp_surface')

    def ok(self):
        global cv
        # Has cb1 a value?
        sel1 = BPSFrame.cb1.getvalue()
        if len(sel1) == 0:
            dialog = Pmw.MessageDialog(BPSFrame,
                title = 'Error',
                defaultbutton = 0,
                buttons = ('OK',),
                message_text = 'ERROR: No ligand selection. Please\n choose one from the dropdownlist.'
            )
            dialog.iconname('Selection Error')
            dialog.bell()
            dialog.activate()
            return 0

        # Has cb2 a value?
        sel2 = BPSFrame.cb2.getvalue()
        if len(sel2) == 0:
            dialog = Pmw.MessageDialog(BPSFrame,
                title = 'Error',
                defaultbutton = 0,
                buttons = ('OK',),
                message_text = 'ERROR: No binding pocket selection. Please\n choose one from the dropdownlist.'
            )
            dialog.iconname('Selection Error')
            dialog.bell()
            dialog.activate()
            return 0

        # Are cb1 and cb2 the same?
        if sel1[0]==sel2[0]:
            dialog = Pmw.MessageDialog(BPSFrame,
                title = 'Error',
                defaultbutton = 0,
                buttons = ('OK',),
                message_text = 'ERROR: Ligand and binding pocket must have different selections'
            )
            dialog.iconname('Selection Error')
            dialog.bell()
            dialog.activate()
            return 0

        # Set orgin and normal vector
        cmd.origin(sel1[0])
        vector = [[  0,  0,  1]]

        # Remap the normal vector to the model space
        d=cmd.get_view(0)
        vector = map( lambda p,v=d: [
            v[0] * p[0] + v[1] * p[1] + v[2]* p[2],
            v[3] * p[0] + v[4] * p[1] + v[5]* p[2],
            v[6] * p[0] + v[7] * p[1] + v[8]* p[2]
            ], vector )

        # Calculate plane coefficients
        # Ax + By + Cz + D = 0
        # D = -Ax -By -Cz
        A = vector[0][0]
        B = vector[0][1]
        C = vector[0][2]
        D = -vector[0][0]*d[12] - vector[0][1]*d[13] - vector[0][2]*d[14]

        atoms = cmd.get_model(sel2[0])

        self.clear_flags(0,1)
        for at in atoms.atom:
            if (A*(at.coord[0])+B*at.coord[1]+C*at.coord[2]+D-BPSFrame.slider.get()>0):
                if not('bp_cutoff' in cmd.get_names("selections")):
                    cmd.select('bp_cutoff','index '+str(at.index))
                else:
                    cmd.select('bp_cutoff','bp_cutoff or index '+str(at.index))
        if ('bp_cutoff' in cmd.get_names("selections")):
            if cv.get()==0:
                cmd.flag('ignore','bp_cutoff','reset')
                cmd.delete('indicate')
                cmd.deselect()
                cmd.show('surface',sel2[0])
            else:
                self.clear_flags(0,0)
                cmd.hide('surface',sel2[0])
                cmd.set('auto_zoom','off')
                cmd.create('bp_surface',sel2[0]+' and not bp_cutoff')
                cmd.set('auto_zoom','on')
                cmd.show('surface','bp_surface')
            cmd.rebuild(sel2[0])

# -*- coding: utf-8 -*-

# -------------------------------------------------------------------------------
# bp_surface.py - Binding pocket surface generator for Pymol version 1.0
# -------------------------------------------------------------------------------

# Copyright (C) 2008 by Lightnir - lightnir@gmail.com

# This script is free software; you can redistribute it and#or modify
# it under the terms of the GNU Library General Public License as
# published by the Free Software Foundation; either version 2 of the
# License, or (at your option) any later version.

# This script is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.

# You should have received a copy of the GNU Library General Public
# License along with this program; if not, write to the
# Free Software Foundation, Inc.,
# 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

from pymol import cmd
from Tkinter import *
import Pmw

BPSFrame = None
selections = []
cv = IntVar()

def __init__(self):
    self.menuBar.addcascademenu('Plugin', 'LightnirsPlugins', 'Lightnirs Plugins',
                                    label='Lightnirs Plugins')
    self.menuBar.addmenuitem('LightnirsPlugins', 'command',
                                 'BPsurface',
                                 label='BP surface',
                                 command = lambda s=self: BPSurface(s))
class BPSurface:
    def __init__(self,app):
        global BPSFrame
        global selections
        global cv
        if BPSFrame==None:
            BPSFrame = Toplevel(width=400, height=300)              # create a root window
            BPSFrame.title("BP surface editor")                     # Set window title
            BPSFrame.resizable(0, 0)                                # Disable window resize
            BPSFrame.geometry('-40+40')                             # Set window placement
            self.get_selections()
            # Create Ligand selection Combobox
            BPSFrame.cb1 = Pmw.ComboBox(BPSFrame,
                label_text = 'Ligand selection:',
                labelpos = 'nw',
                scrolledlist_items = selections
            )
            BPSFrame.cb1.grid(row=0, column=0,columnspan=2)

            # Create binding pocket Combobox
            BPSFrame.cb2 = Pmw.ComboBox(BPSFrame,
                label_text = 'Binding pocket selection:',
                labelpos = 'nw',
                scrolledlist_items = selections
            )
            BPSFrame.cb2.grid(row=1, column=0,columnspan=2)

            # Create Slider
            BPSFrame.slider = Scale(BPSFrame, from_=-20, to=20, tickinterval=10, resolution=0.5, label='Z axis cutoff', orient=HORIZONTAL)
            BPSFrame.slider.grid(row=2, column=0, columnspan=2, sticky=N+W+S+E,)

            # Add Checkbox
            cbx = Checkbutton(BPSFrame, text="Create new object", variable=cv, onvalue="1", offvalue="0")
            cbx.select()
            cbx.grid(row=3,column=0, columnspan=2, sticky=W)

            # Add some buttons
            btn1=Button(BPSFrame, text='Refresh', padx=0, pady=0, command = self.refresh_list)
            btn1.grid(row=4, column=0, sticky=N+W+S+E, padx=0, pady=0,ipadx=0,ipady=0)
            btn2=Button(BPSFrame, text='Clear', padx=0, pady=0, command = lambda r=1,d=1: self.clear_flags(r,d))
            btn2.grid(row=4, column=1, sticky=N+W+S+E, padx=0, pady=0,ipadx=0,ipady=0)
            btn3=Button(BPSFrame, text='Apply', padx=0, pady=0, command = self.ok)
            btn3.grid(row=5, column=0, sticky=N+W+S+E,)
            btn4=Button(BPSFrame, text='Close', padx=0, pady=0, command = self.myHide)
            btn4.grid(row=5, column=1, sticky=N+W+S+E)

            # create callback to prevent window kill
            BPSFrame.protocol("WM_DELETE_WINDOW", self.myHide)
            BPSFrame.mainloop()
        else:
            if BPSFrame.state() == "normal":
                self.myHide()
            elif BPSFrame.state() == "withdrawn":
                self.myShow()

    def myShow(self):
        self.refresh_list()
        BPSFrame.deiconify()

    def myHide(self):
        BPSFrame.withdraw()

    def get_selections(self):
        global selections
        selections = []
        for item in cmd.get_names("selections"):
            if item[0]<>"_":
                selections.append(item)

    def refresh_list(self):
        global selections

        # Save old values
        old = selections
        try:
            old_cb1=BPSFrame.cb1.getvalue()[0]
        except:
            old_cb1=''
        try:
            old_cb2=BPSFrame.cb2.getvalue()[0]
        except:
            old_cb2=''
        self.get_selections()

        # Update selectionlist if needed
        if not(old==selections):
            BPSFrame.cb1.setlist(selections)
            BPSFrame.cb2.setlist(selections)

        # Is ligand selection still available?
        if not(old_cb1 in selections):
            # No - clear the combobox entry
            BPSFrame.cb1.component('entryfield').clear()
        else:
            # Yes - set the old value
            BPSFrame.cb1.selectitem(old_cb1,setentry=1)

        # Is binding pocket selection still available?
        if not(old_cb2 in selections):
            # No - clear the combobox entry
            BPSFrame.cb2.component('entryfield').clear()
        else:
            # Yes - set the old value
            BPSFrame.cb2.selectitem(old_cb2,setentry=1)

    def clear_flags(self,r,d):
        cmd.flag('ignore','all','clear')
        if r==1:
            cmd.rebuild('all')
        if d==1:
            cmd.delete('bp_cutoff')
            cmd.delete('bp_surface')

    def ok(self):
        global cv
        # Has cb1 a value?
        sel1 = BPSFrame.cb1.getvalue()
        if len(sel1) == 0:
            dialog = Pmw.MessageDialog(BPSFrame,
                title = 'Error',
                defaultbutton = 0,
                buttons = ('OK',),
                message_text = 'ERROR: No ligand selection. Please\n choose one from the dropdownlist.'
            )
            dialog.iconname('Selection Error')
            dialog.bell()
            dialog.activate()
            return 0

        # Has cb2 a value?
        sel2 = BPSFrame.cb2.getvalue()
        if len(sel2) == 0:
            dialog = Pmw.MessageDialog(BPSFrame,
                title = 'Error',
                defaultbutton = 0,
                buttons = ('OK',),
                message_text = 'ERROR: No binding pocket selection. Please\n choose one from the dropdownlist.'
            )
            dialog.iconname('Selection Error')
            dialog.bell()
            dialog.activate()
            return 0

        # Are cb1 and cb2 the same?
        if sel1[0]==sel2[0]:
            dialog = Pmw.MessageDialog(BPSFrame,
                title = 'Error',
                defaultbutton = 0,
                buttons = ('OK',),
                message_text = 'ERROR: Ligand and binding pocket must have different selections'
            )
            dialog.iconname('Selection Error')
            dialog.bell()
            dialog.activate()
            return 0

        # Set orgin and normal vector
        cmd.origin(sel1[0])
        vector = [[  0,  0,  1]]

        # Remap the normal vector to the model space
        d=cmd.get_view(0)
        vector = map( lambda p,v=d: [
            v[0] * p[0] + v[1] * p[1] + v[2]* p[2],
            v[3] * p[0] + v[4] * p[1] + v[5]* p[2],
            v[6] * p[0] + v[7] * p[1] + v[8]* p[2]
            ], vector )

        # Calculate plane coefficients
        # Ax + By + Cz + D = 0
        # D = -Ax -By -Cz
        A = vector[0][0]
        B = vector[0][1]
        C = vector[0][2]
        D = -vector[0][0]*d[12] - vector[0][1]*d[13] - vector[0][2]*d[14]

        atoms = cmd.get_model(sel2[0])

        self.clear_flags(0,1)
        for at in atoms.atom:
            if (A*(at.coord[0])+B*at.coord[1]+C*at.coord[2]+D-BPSFrame.slider.get()>0):
                if not('bp_cutoff' in cmd.get_names("selections")):
                    cmd.select('bp_cutoff','index '+str(at.index))
                else:
                    cmd.select('bp_cutoff','bp_cutoff or index '+str(at.index))
        if ('bp_cutoff' in cmd.get_names("selections")):
            if cv.get()==0:
                cmd.flag('ignore','bp_cutoff','reset')
                cmd.delete('indicate')
                cmd.deselect()
                cmd.show('surface',sel2[0])
            else:
                self.clear_flags(0,0)
                cmd.hide('surface',sel2[0])
                cmd.set('auto_zoom','off')
                cmd.create('bp_surface',sel2[0]+' and not bp_cutoff')
                cmd.set('auto_zoom','on')
                cmd.show('surface','bp_surface')
            cmd.rebuild(sel2[0])

Put the script bp_surface.py in the pymol startup directory, that is under Windows:

C:\Program Files\DeLano Scientific\PyMOL\modules\pmg_tk\startup
or under Linux (as in my case):

/usr/share/pymol/modules/pmg_tk/startup
Now lets get a model to work with. I'm in a right mood for some complexed dopamine so lets use the 2a3r PDB entry. The plugin requires two selections to work. First the ligand and second - the surrounding residues. In this enzyme model there are two L-dopamine molecule. Lets select one of them by typing:

select ligand, resn LDP and resi 201
Now for the second selection.
select pocket, byres(ligand around 8)
After this we are ready to work with the plugin. After starting it (Plugins -> Lightnirs plugins -> BP surface) we need to select the right selections in the comboboxes. The slider is for adjusting the cutting plane on the Z axis. With it we can make the ligand look more or less sunken in the surrounding residue surface. For example, by setting positive values we move the cutting plane towards the viewer and the ligand looks more sunken. The Refresh button updates the selections-list that is available in the two comboboxes (just in case you made a new selection after starting the plugin). The Clear button removes all modifications made by this plugin. If the "Create new object" checkbox is checked(default) then there will be created a new object for the surface and the surface will be closed, otherwise you'll get a partial surface.

Here's how it looks for the L-dopamine (after changing some color and representation settings):

Notice that this plugin creates a "bp_cutoff" selection which is a subset of the "apple flesh" and can be used to quickly hide the residues that are in the way.

Here's the view from the side. Notice the clean cut - it's only one part of the binding site.

And here is the second one.


Showing a binding site is now as easy as cutting an apple. I'm thinking about putting this script on the Pymol wiki, but first I'll do some more testing on it. Any suggestions on this are welcome ]:)

UPDATE: There's a new version of my plugin here