Hi PyMol-Community,
I have been looking for a PyMOL command to calculate the ribose sugar pucker
information (in DNA) but have not found anything. Thus, I began writing a
simple Python script that is supposed to take a selection, determine whether or
not it contains a ribose sugar ring, and then calculate the phase, amplitude,
and pucker orientation for each nucleotide residue within the selection. The
completed script works just fine (see the script below) but it slows down
significantly in performance when calling cmd.get_dihedral (at least, my crude
debugging method using print statements has identified the culprit to be the
get_dihedral command) with a LARGE protein-DNA complex (160,000 atoms) loaded.
However, if I only load the DNA (with no protein), the performance is lightning
fast. Can anybody help me improve the overall performance or suggest an
alternative? Thanks in advance!
Sean
=========================================================
from pymol.cgo import * # get constants
from math import *
from pymol import cmd
def pucker(selection):
# Author: Sean Law
# Institute: Michigan State University
# E-mail: s...@msu.edu
obj_name = selection
objects = cmd.get_names(selection)
for obj in objects:
if (obj == selection):
obj_name=obj
sel=cmd.get_model(selection)
first=1
old=" "
oldchain=" "
residue={}
count=0
for atom in sel.atom:
new=atom.chain+" "+str(atom.resi)
newchain=atom.chain+" "+atom.segi
if (not (oldchain == newchain) and first):
print " " #Blank line to separate chain output
print "%6s %6s %8s Residue" % ("Phase", "Amp", "Pucker")
if (not(new == old) and (not first)):
#Check that all 5 atoms exist
if(len(residue) == 5):
#Calculate pucker
info = pseudo(residue)
print info+" "+old
else:
print "There is no sugar in this residue"
if (not (oldchain == newchain)):
print " " #Blank line to separate chain output
print "%6s %6s %8s Residue" % ("Phase", "Amp", "Pucker")
#Clear values
residue={}
#Store new value
store_atom(atom,residue,obj_name)
else:
store_atom(atom,residue,obj_name)
first=0
old=new
oldchain=newchain
#Final Residue
#Calculate dihedrals for final residue
if (len(residue) == 5):
#Calculate pucker for final residue
info = pseudo(residue)
print info+" "+old
else:
print "There is no sugar in this residue"
return
def sele_exists(sele):
return sele in cmd.get_names("selections");
def pseudo(residue):
other=2*(sin(math.radians(36.0))+sin(math.radians(72.0)))
theta0=cmd.get_dihedral(residue['C4*'],residue['O4*'],residue['C1*'],residue['C2*'])
theta1=cmd.get_dihedral(residue['O4*'],residue['C1*'],residue['C2*'],residue['C3*'])
theta2=cmd.get_dihedral(residue['C1*'],residue['C2*'],residue['C3*'],residue['C4*'])
theta3=cmd.get_dihedral(residue['C2*'],residue['C3*'],residue['C4*'],residue['O4*'])
theta4=cmd.get_dihedral(residue['C3*'],residue['C4*'],residue['O4*'],residue['C1*'])
#phase=atan2((theta4+theta1)-(theta3+theta0),2*theta2*(sin(math.radians(36.0))+sin(math.radians(72.0))))
phase=atan2((theta4+theta1)-(theta3+theta0),theta2*other)
amplitude=theta2/cos(phase)
phase=math.degrees(phase)
if (phase < 0):
phase=phase+360 # 0 <= Phase < 360
#Determine pucker
if (phase < 36):
pucker = "C3'-endo"
elif (phase < 72):
pucker = "C4'-exo"
elif (phase <108):
pucker = "O4'-endo"
elif (phase < 144):
pucker = "C1'-exo"
elif (phase < 180):
pucker = "C2'-endo"
elif (phase < 216):
pucker = "C3'-exo"
elif (phase < 252):
pucker = "C4'-endo"
elif (phase < 288):
pucker = "O4'-exo"
elif (phase < 324):
pucker = "C1'-endo"
elif (phase < 360):
pucker = "C2'-exo"
else:
pucker = "Phase is out of range"
info = "%6.2f %6.2f %8s" % (phase, amplitude, pucker)
return info
def store_atom(atom,residue,obj_name):
if (atom.name == "O4'" or atom.name == "O4*"):
residue['O4*'] = str(atom_sele(atom,obj_name))
elif (atom.name == "C1'" or atom.name == "C1*"):
residue['C1*'] = str(atom_sele(atom,obj_name))
elif (atom.name == "C2'" or atom.name == "C2*"):
residue['C2*'] = str(atom_sele(atom,obj_name))
elif (atom.name == "C3'" or atom.name == "C3*"):
residue['C3*'] = str(atom_sele(atom,obj_name))
elif (atom.name == "C4'" or atom.name == "C4*"):
residue['C4*'] = str(atom_sele(atom,obj_name))
return
def atom_sele(atom,obj_name):
atom_sele = ""
if (not (obj_name == "")):
atom_sele = atom_sele+obj_name+" & "
if (not (atom.segi == "")):
atom_sele = atom_sele+"segi "+atom.segi+" & "
if (not (atom.chain == "")):
atom_sele = atom_sele+"chain "+atom.chain+" & "
if (not (atom.resn == "")):
atom_sele = atom_sele+"resn "+atom.resn+" & "
if (not (atom.resi == "")):
atom_sele = atom_sele+"resi "+str(atom.resi)+" & "
if (not (atom.name == "")):
atom_sele = atom_sele+"name "+atom.name
return atom_sele
cmd.extend("pucker",pucker)
_________________________________________________________________
Reinvent how you stay in touch with the new Windows Live Messenger.
http://go.microsoft.com/?linkid=9650731
