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)
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.
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.
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 def 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