Rotamer Toggle

DESCRIPTION

Backbone-Dependent Rotamer library (Dunbrack, Cohen ; see ref) is imported into pymol giving access to this information. There are a number of different ways to use the data, I've only implemented a few as well as added extra functions that seemed useful.

• colorRotamers - color rotamers by closest matching rotamer angles from database; i.e. color by how common each rotamer of selection, blue - red (least to most common).
• Rotamer Menu - an added menu into menu.py, which displays the most common rotamers for the given(clicked) residue; you can also set the residue any of the common rotamers as well
• set_rotamer - routine called by above menu, but can be called manually to set a specific residues side-chain angles
• set_phipsi - set all phi,psi angles of given selection to given angles (useful for creating secondary structures)

SETUP

To setup a rotamer menu inside the residue menu (default windows pymol installation):

• copy rotamers.py to C:/Program Files/DeLano Scientific/PyMol/modules/pymol/rotamers.py
• copy mymenu.py to C:/Program Files/DeLano Scientific/PyMol/modules/pymol/menu.py (WARNING : overwrites default menu.py - use at your own risk)
• copy bbdep02.May.sortlib to C:/Program Files/DeLano Scientific/PyMol/modules/pymol/bbdep02.May.sortlib (or newer version of sorted bbdep)

This is only one possible way to do this, I am sure there are many others. I'm not going to post the bbdep, but there is a link in the References section to Dunbrack's download page (get the "sorted" lib)

Of course you can simply import rotamers.py and use the functions manually.

NOTES / STATUS

This should not have any OS specific stuff in it yet, but has only been tested on Windows. I'll test it on linux in the next few days. The way it's setup now, when you import rotamers , it will automatically read-in the rotamer database; this may not be what you want.

USAGE

colorRotamers selection
set_rotamer selection, chi1_angle [,chi2_angle] [,chi3_angle] [,chi4_angle]
set_phipsi selection phi_angle, psi_angle

EXAMPLES

  colorRotamers chain A
  set_rotamer resi 40, -60,-40   (only set chi1,chi2 angles)
  set_phipsi resi 10-40, -60,-60 (create an alpha-helical-like section)
  

REFERENCES

Dunbrack and Cohen. Protein Science 1997

Dunbrack Lab Page (Contains backbone-dependent library)

Scripts (Rotamers.py ; MyMenu.py)

Rotamers.py

##################################################################
# File:          Rotamers.py
# Author:        Dan Kulp
# Creation Date: 6/8/05
#
# Notes:
#     Incorporation of Rotamer library
#     readRotLib() - fills rotdat; 
#        indexed by "RES:PHI_BIN:PSI_BIN".
#
#     Three main functions:
#     1. colorRotamers - colors according
#          to rotamer probablitity
#     2. getBins(sel)
#           phi,psi bin for rotamer
#     3. set_rotamer - set a side-chain 
#           to a specific rotamer	
#
#     To setup a rotamer menu in the 
#   right click, under "Residue"
#        1. cp mymenu.py modules/pymol/menu.py
#        2. cp rotamers.py modules/pymol/rotamers.py (update ROTLIB)
#
# Requirements:
#  set ROTLIB to path for rotamer library
# Reference: 
#  Dunbrack and Cohen. Protein Science 1997
####################################################################

import colorsys,sys
import string
import re
import editing
import os
import cmd
import math

# Path for library
ROTLIB=os.environ['PYMOL_PATH']+"/modules/pymol/bbdep02.May.sortlib"

# Place for library in memory..
rotdat = {}

def readRotLib():	
    # Column indexes in rotamer library..
    RES  = 0
    PHI  = 1
    PSI  = 2
    PROB = 8
    CHI1 = 9
    CHI2 = 10
    CHI3 = 11
    CHI4 = 12

    if os.path.exists(ROTLIB):
		print "File exists: "+ROTLIB
		input = open(ROTLIB, 'r')
		for line in input:

	  	    # Parse by whitespace (I believe format is white space and not fixed-width columns)
		    dat = re.split("\s+",line)

		    # Add to rotamer library in memory : 
		    #   key format       RES:PHI_BIN:PSI_BIN
		    #   value format     PROB, CHI1, CHI2, CHI3, CHI4
		    try:
		        rotdat[dat[RES]+":"+dat[PHI]+":"+dat[PSI]].append([ dat[PROB], dat[CHI1], dat[CHI2], dat[CHI3], dat[CHI4] ])
		    except KeyError:
			rotdat[dat[RES]+":"+dat[PHI]+":"+dat[PSI]] = [ [ dat[PROB], dat[CHI1], dat[CHI2], dat[CHI3], dat[CHI4] ] ]

		    
    else:
	print "Couldn't find Rotamer library"


# Atoms for each side-chain angle for each residue
CHIS = {}
CHIS["ARG"] = [ ["N","CA","CB","CG" ],
		["CA","CB","CG","CD" ],
		["CB","CG","CD","NE" ],
		["CG","CD","NE","CZ" ]
	      ]

CHIS["ASN"] = [ ["N","CA","CB","CG" ],
		["CA","CB","CG","OD2" ]
	      ]

CHIS["ASP"] = [ ["N","CA","CB","CG" ],
		["CA","CB","CG","OD1" ]
	      ]
CHIS["CYS"] = [ ["N","CA","CB","CG" ]
	      ]
CHIS["GLN"] = [ ["N","CA","CB","CG" ],
		["CA","CB","CG","CD" ],
		["CB","CG","CD","OE1"]
	      ]

CHIS["GLU"] = [ ["N","CA","CB","CG" ],
		["CA","CB","CG","CD" ],
		["CB","CG","CD","OE1"]
	      ]

CHIS["HIS"] = [ ["N","CA","CB","CG" ],
		["CA","CB","CG","ND1"]
	      ]

CHIS["ILE"] = [ ["N","CA","CB","CG1" ],
		["CA","CB","CG1","CD1" ]
	      ]

CHIS["LEU"] = [ ["N","CA","CB","CG" ],
		["CA","CB","CG","CD1" ]
	      ]

CHIS["LYS"] = [ ["N","CA","CB","CG" ],
		["CA","CB","CG","CD" ],
		["CB","CG","CD","CE"],
		["CG","CD","CE","NZ"]
	      ]

CHIS["MET"] = [ ["N","CA","CB","CG" ],
		["CA","CB","CG","SD" ],
		["CB","CG","SD","CE"]
	      ]

CHIS["PHE"] = [ ["N","CA","CB","CG" ],
		["CA","CB","CG","CD1" ]
	      ]

CHIS["PRO"] = [ ["N","CA","CB","CG" ],
		["CA","CB","CG","CD" ]
	      ]

CHIS["SER"] = [ ["N","CA","CB","OG" ]
	      ]

CHIS["THR"] = [ ["N","CA","CB","OG1" ]
	      ]
CHIS["TYR"] = [ ["N","CA","CB","CG" ],
		["CA","CB","CG","CD1" ]
	      ]

CHIS["VAL"] = [ ["N","CA","CB","CG1" ]
	      ]

# Color Rotamer by side-chain angle position
#  'bin' side-chain angles into closest
def colorRotamers(sel):
    doRotamers(sel)

# Utility function, to set phi,psi angles for a given selection
# Note: Cartoon, Ribbon functionality will not display correctly after this
def set_phipsi(sel, phi,psi):
    doRotamers(sel,angles=[phi,psi],type="set")

# Set a rotamer, based on a selection, a restype and chi angles
def set_rotamer(sel, chi1, chi2=0,chi3=0,chi4=0):
    at = cmd.get_model("byres ("+sel+")").atom[0]

    list = [chi1,chi2,chi3,chi4]
    for i in range(0,len(CHIS[at.resn])):
	print "Setting Chi"+str(i+1)+" to "+str(list[i])
        editing.set_dihedral(sel + ' and name '+CHIS[at.resn][i][0],		    
		             sel + ' and name '+CHIS[at.resn][i][1],		    
		             sel + ' and name '+CHIS[at.resn][i][2],		    
		             sel + ' and name '+CHIS[at.resn][i][3], list[i])		    

    # Remove some objects that got created
    cmd.delete("pk1")
    cmd.delete("pk2")
    cmd.delete("pkmol")

# Get Phi,Psi bins for given selection
# WARNING:  assume selection is single residue (will only return first residue bins)
def getBins(sel):
    return doRotamers(sel, type="bins")

# Specific comparison operator for rotamer prob data 
def mycmp(first, second):
	return cmp( first[1], second[1])

# Color Ramp...
def rot_color(vals): 
	nbins = 10
	vals.sort(mycmp)
	print "End sort: "+str(len(vals))+" : "+str(nbins)


	# Coloring scheme...
	i = 0
	j = 0
	rgb = [0.0,0.0,0.0]
	sel_str = ""
	while i < len(vals):
		if int(len(vals)/nbins) == 0 or i % int(len(vals)/nbins) == 0:
		      hsv = (colorsys.TWO_THIRD - colorsys.TWO_THIRD * float(j) / (nbins-1), 1.0, 1.0)

		      #convert to rgb and append to color list
		      rgb = colorsys.hsv_to_rgb(hsv[0],hsv[1],hsv[2])
		      if j < nbins-1:
		              j += 1	

		cmd.set_color("RotProbColor"+str(i), rgb)
		cmd.color("RotProbColor"+str(i), str(vals[i][0]))
		i += 1


# Main function                           
def doRotamers(sel,angles=[], type="color"):                           

	# Read in Rotamer library if not already done
	if len(rotdat) == 0:
		readRotLib()

	# Set up some variables..
	residues = ['dummy']  # Keep track of residues already done
	probs = []            # probability of each residue conformation
	phi = 0               # phi,psi angles of current residue
	psi = 0

	# Get atoms from selection
	atoms = cmd.get_model("byres ("+sel+")")

        # Loop through atoms in selection		
	for at in atoms.atom:
	    try:
	       # Don't process Glycines or Alanines
	       if not (at.resn == 'GLY' or at.resn == 'ALA'):
	        if not at.chain+":"+at.resn+":"+at.resi in residues:
	            residues.append(at.chain+":"+at.resn+":"+at.resi)

		    # Check for a null chain id (some PDBs contain this) 
		    unit_select = ""
		    if not at.chain == "":
			unit_select = "chain "+str(at.chain)+" and "

		    # Define selections for residue i-1, i and i+1    
		    residue_def = unit_select+'resi '+str(at.resi)
  		    residue_def_prev = unit_select+'resi '+str(int(at.resi)-1)
		    residue_def_next = unit_select+'resi '+str(int(at.resi)+1)

	            # Compute phi/psi angle
		    phi = cmd.get_dihedral(residue_def+' and name CB',residue_def+' and name CA',residue_def+' and name N',residue_def_prev+' and name C')
		    psi = cmd.get_dihedral(residue_def+' and name O',residue_def+' and name C',residue_def+' and name CA',residue_def+' and name CB')
		    if type == "set":
			    print "Changing "+at.resn+str(at.resi)+" from "+str(phi)+","+str(psi)+" to "+str(angles[0])+","+str(angles[1])
			    cmd.set_dihedral(residue_def+' and name CB',residue_def+' and name CA',residue_def+' and name N',residue_def_prev+' and name C',angles[0])
			    cmd.set_dihedral(residue_def+' and name O',residue_def+' and name C',residue_def+' and name CA',residue_def+' and name CB', angles[1])
			    continue
						
		    # Find correct 10x10 degree bin		        		
		    phi_digit = abs(int(phi)) - abs(int(phi/10)*10)
		    psi_digit = abs(int(psi)) - abs(int(psi/10)*10)
		    
		    # Remember sign of phi,psi angles
	    	    phi_sign = 1
		    if phi < 0:    phi_sign = -1

		    psi_sign = 1
		    if psi < 0:    psi_sign = -1

		    # Compute phi,psi bins
  		    phi_bin = int(math.floor(abs(phi/10))*10*phi_sign)
		    if phi_digit >= 5:    phi_bin = int(math.ceil(abs(phi/10))*10*phi_sign)

		    psi_bin = int(math.floor(abs(psi/10))*10*psi_sign)
		    if psi_digit >= 5:    psi_bin = int(math.ceil(abs(psi/10))*10*psi_sign)

	            print "Given "+at.resn+":"+at.resi+" PHI,PSI ("+str(phi)+","+str(psi)+") : bin ("+str(phi_bin)+","+str(psi_bin)+")"

		
		    # Get current chi angle measurements
		    chi = []
		    for i in range(0,len(CHIS[at.resn])):
		       chi.append(cmd.get_dihedral(residue_def + ' and name '+CHIS[at.resn][i][0],		    
		  			             residue_def + ' and name '+CHIS[at.resn][i][1],		    
					             residue_def + ' and name '+CHIS[at.resn][i][2],		    
					             residue_def + ' and name '+CHIS[at.resn][i][3]))		    
		    print "CHIs: "+str(chi)
		    if type == 'bins':
		         return [at.resn, phi_bin,psi_bin]

		    # Compute probabilities for given chi angles
                    prob = 0
		    prob_box = 22 		    
		    for item in range(0,len(rotdat[at.resn+":"+str(phi_bin)+":"+str(psi_bin)])):
			print "Rotamer from db: "+str(rotdat[at.resn+":"+str(phi_bin)+":"+str(psi_bin)][item])
			if chi[0]:
			    if chi[0] >= float(rotdat[at.resn+":"+str(phi_bin)+":"+str(psi_bin)][item][1]) - (prob_box/2) and \
				chi[0] <= float(rotdat[at.resn+":"+str(phi_bin)+":"+str(psi_bin)][item][1]) + (prob_box/2):
				if len(chi) == 1:
					prob = rotdat[at.resn+":"+str(phi_bin)+":"+str(psi_bin)][item][0]
					break
				if chi[1] >= float(rotdat[at.resn+":"+str(phi_bin)+":"+str(psi_bin)][item][2]) - (prob_box/2) and \
				 float(chi[1] <= rotdat[at.resn+":"+str(phi_bin)+":"+str(psi_bin)][item][2]) + (prob_box/2):
					if len(chi) == 2:
					    prob = rotdat[at.resn+":"+str(phi_bin)+":"+str(psi_bin)][item][0]
					    break
					if chi[2] >= float(rotdat[at.resn+":"+str(phi_bin)+":"+str(psi_bin)][item][3]) - (prob_box/2) and \
					   float(chi[2] <= rotdat[at.resn+":"+str(phi_bin)+":"+str(psi_bin)][item][3]) + (prob_box/2):
    					    if len(chi) == 3:
					        prob = rotdat[at.resn+":"+str(phi_bin)+":"+str(psi_bin)][item][0]
					        break
					    if chi[3] >= float(rotdat[at.resn+":"+str(phi_bin)+":"+str(psi_bin)][item][4]) - (prob_box/2) and \
					       float(chi[3] <= rotdat[at.resn+":"+str(phi_bin)+":"+str(psi_bin)][item][4]) + (prob_box/2):
					        prob = rotdat[at.resn+":"+str(phi_bin)+":"+str(psi_bin)][item][0]
					        break

			
		    print "PROB OF ROTAMER: "+str(prob)
		    print "---------------------------"
		    probs.append([residue_def, prob])

	    except:
#		probs.append([residue_def, -1])
		print "Exception found"
		continue

	# Color according to rotamer probability
	rot_color(probs)	    

	




readRotLib()

cmd.extend('set_phipsi',set_phipsi)
cmd.extend('set_rotamer',set_rotamer)
cmd.extend('colorRotamers',colorRotamers)

MyMenu.py

Since menu.py is copyrighted I can't post my edited version, but you can create it very simply by adding these two peices of code

1. In the "pick_option(title,s,object=0)" function of menu.py add the following code after the first "result =" statement

# Edit dwkulp 6/11/05 , add a rotamer menu to residue menu
   if title == 'Residue':
	result.extend([[ 1, 'rotamers'    , rotamer_menu(s)]])

2. At the end of the file add this:

###############################################
# Dan Kulp
# Added Rotamer list to residue menu..
# rotamer.py must be importable (I placed it in 
# the same directory as menu.py)
###############################################

import rotamers


def rotamer_menu(s):
	# Check for rotamer library being loaded
	if not rotamers.rotdat:
	     return [ [2, "Must run rotamers.py first",'']]

	# Check for valid rotamer residue..
	res = cmd.get_model("byres ("+s+")").atom[0].resn
        if not res in rotamers.CHIS.keys():
	    return [ [2, "Residue: "+res+" not known sidechain or does not have rotamers", '']]

	# Get PHI,PSI bins for rotamer (also prints out current phi,psi, chi1,chi2,chi3,chi4)
	bins = rotamers.doRotamers(s,type='bins')

	# Add a title to the menu
	result = [ [2, bins[0]+' Rotamers in bin ('+str(bins[1])+','+str(bins[2])+')','' ], [1, ':::PROB,CHI1,CHI2,CHI3,CHI4:::','']]

        # Grab the entries for this residue and phi,psi bins
	match_rotamers = rotamers.rotdat[bins[0]+":"+str(bins[1])+":"+str(bins[2])]

	# Set max number of rotamers to display (probably should be somewhere 'higher up' in the code)
	max_rotamers = 10


	if len(match_rotamers) < max_rotamers:
	    max_rotamers = len(match_rotamers)

	# Create menu entry for each possible rotamer
        for item in range(0,max_rotamers):
             result.append( [ 1, str(match_rotamers[item]), 'rotamers.set_rotamer("'+s+'","'\
										    +str(match_rotamers[item][1])+'","'\
										    +str(match_rotamers[item][2])+'","'\
										    +str(match_rotamers[item][3])+'","'\
										    +str(match_rotamers[item][4])+'")'])
	return result