This is a read-only mirror of pymolwiki.org

Difference between revisions of "Pucker"

From PyMOL Wiki
Jump to navigation Jump to search
m (17 revisions)
 
(10 intermediate revisions by 2 users not shown)
Line 1: Line 1:
 +
{{Infobox script-repo
 +
|type      = script
 +
|filename  = pucker.py
 +
|author    = [[User:Slaw|Sean M. Law]]
 +
|license  = -
 +
}}
 +
 
===DESCRIPTION===  
 
===DESCRIPTION===  
 
'''pucker.py''' is a PyMol script that returns the sugar pucker information (phase, amplitude, pucker) for a given selection.
 
'''pucker.py''' is a PyMol script that returns the sugar pucker information (phase, amplitude, pucker) for a given selection.
  
−
This script does NOT uses its own dihedral calculation scheme rather than the get_dihedral command.  Thus, it is lightning fast!
+
This script uses its own dihedral calculation scheme rather than the get_dihedral command.  Thus, it is lightning fast!
  
 
If a selection does not contain any ribose sugars then an error message is returned.   
 
If a selection does not contain any ribose sugars then an error message is returned.   
Line 8: Line 15:
 
===USAGE===
 
===USAGE===
  
−
load the script using the "run" command
+
load the script using the [[run]] command
  
 +
<source lang="python">
 
pucker selection
 
pucker selection
 +
#Calculates the sugar information for the given selection
 +
 +
pucker selection, state=1
 +
#Calculates the sugar information for the given selection in state 1
 +
 +
pucker selection, state=0
 +
#Iterates over all states and calculates the sugar information for the given selection
 +
</source>
  
 
===EXAMPLES===
 
===EXAMPLES===
 +
<source lang="python">
 +
#fetch 1BNA.pdb
 +
fetch 1bna
  
−
pucker all
+
#Select DNA only
 +
#Otherwise, you will get an error for water not having sugars
 +
select DNA, not solvent
 +
 
 +
#Execute pucker command
 +
pucker DNA
 +
 
 +
#The output should look like this
 +
Phase    Amp    Pucker  Residue
 +
161.22  55.32  C2'-endo  A 1
 +
139.52  41.67  C1'-exo  A 2
 +
92.82  38.31  O4'-endo  A 3
 +
166.35  48.47  C2'-endo  A 4
 +
128.57  46.30  C1'-exo  A 5
 +
126.92  49.75  C1'-exo  A 6
 +
101.30  47.32  O4'-endo  A 7
 +
115.62  49.23  C1'-exo  A 8
 +
140.44  46.37  C1'-exo  A 9
 +
145.79  53.36  C2'-endo  A 10
 +
147.47  47.04  C2'-endo  A 11
 +
113.80  51.69  C1'-exo  A 12
 +
 +
Phase    Amp    Pucker  Residue
 +
153.24  43.15  C2'-endo  B 13
 +
128.49  45.01  C1'-exo  B 14
 +
67.74  43.84  C4'-exo  B 15
 +
149.33  41.01  C2'-endo  B 16
 +
169.27  42.30  C2'-endo  B 17
 +
147.03  42.30  C2'-endo  B 18
 +
116.29  47.52  C1'-exo  B 19
 +
129.62  49.92  C1'-exo  B 20
 +
113.61  42.93  C1'-exo  B 21
 +
156.34  50.98  C2'-endo  B 22
 +
116.89  44.21  C1'-exo  B 23
 +
34.70  45.55  C3'-endo  B 24
 +
</source>
  
  
Line 24: Line 78:
 
from pymol import cmd
 
from pymol import cmd
  
−
def pucker(selection):
+
def pucker(selection,state=1):
  
 
# Author: Sean Law
 
# Author: Sean Law
−
# Institute: Michigan State University
+
# Institute: University of Michigan
−
# E-mail: slaw@msu.edu
+
# E-mail: seanlaw@umich.edu
−
+
 
−
sel=cmd.get_model(selection)
+
#Comparison to output from 3DNA is identical
 +
#Note that the 3DNA output has the second chain
 +
#printed out in reverse order
 +
state=int(state)
 +
if state == 0:
 +
for state in range(1,cmd.count_states()+1):
 +
sel=cmd.get_model(selection,state)
 +
if state == 1:
 +
print " " #Blank line to separate chain output
 +
print "%5s  %6s  %6s  %8s  Residue" % ("State","Phase","Amp", "Pucker")
 +
get_pucker(sel,iterate=state)
 +
else:
 +
sel=cmd.get_model(selection,state)
 +
get_pucker(sel,iterate=0)
 +
return
 +
 
 +
def get_pucker (sel,iterate=0):
 +
iterate=int(iterate)
 
first=1
 
first=1
 
old=" "
 
old=" "
Line 40: Line 111:
 
new=atom.chain+" "+str(atom.resi)
 
new=atom.chain+" "+str(atom.resi)
 
newchain=atom.chain+" "+atom.segi
 
newchain=atom.chain+" "+atom.segi
−
if (not (oldchain == newchain) and first):
+
if oldchain != newchain and first:
−
print " " #Blank line to separate chain output
+
if iterate == 0:
−
print "%6s  %6s  %8s  Residue" % ("Phase", "Amp", "Pucker")
+
print " " #Blank line to separate chain output
−
if (not(new == old) and (not first)):
+
print "%6s  %6s  %8s  Residue" % ("Phase", "Amp", "Pucker")
 +
if new != old and not first:
 
#Check that all 5 atoms exist
 
#Check that all 5 atoms exist
−
if(len(residue) == 15):
+
if len(residue) == 15:
 
#Construct planes
 
#Construct planes
 
get_dihedrals(residue,theta)
 
get_dihedrals(residue,theta)
Line 53: Line 125:
 
else:
 
else:
 
print "There is no sugar in this residue"
 
print "There is no sugar in this residue"
−
if (not (oldchain == newchain)):
+
if oldchain != newchain:
 
print " " #Blank line to separate chain output
 
print " " #Blank line to separate chain output
 
print "%6s  %6s  %8s  Residue" % ("Phase", "Amp", "Pucker")
 
print "%6s  %6s  %8s  Residue" % ("Phase", "Amp", "Pucker")
Line 69: Line 141:
 
#Final Residue
 
#Final Residue
 
#Calculate dihedrals for final residue
 
#Calculate dihedrals for final residue
−
if (len(residue) == 15):
+
if len(residue) == 15:
 
#Construct planes
 
#Construct planes
 
get_dihedrals(residue,theta)
 
get_dihedrals(residue,theta)
 
#Calculate pucker for final residue
 
#Calculate pucker for final residue
 
info = pseudo(residue,theta)
 
info = pseudo(residue,theta)
−
print info+"  "+old
+
if iterate == 0:
 +
print info+"  "+old
 +
else:
 +
print "%5d  %s" % (iterate,info+"  "+old)
 
else:
 
else:
 
print "There is no sugar in this residue"
 
print "There is no sugar in this residue"
Line 89: Line 164:
 
amplitude=theta['2']/cos(phase)
 
amplitude=theta['2']/cos(phase)
 
phase=math.degrees(phase)
 
phase=math.degrees(phase)
−
if (phase < 0):
+
if phase < 0:
−
phase=phase+360 # 0 <= Phase < 360
+
phase+=360 # 0 <= Phase < 360
 
#Determine pucker
 
#Determine pucker
−
if (phase < 36):
+
if phase < 36:
 
pucker = "C3'-endo"
 
pucker = "C3'-endo"
−
elif (phase < 72):
+
elif phase < 72:
 
pucker = "C4'-exo"
 
pucker = "C4'-exo"
−
elif (phase <108):
+
elif phase < 108:
 
pucker = "O4'-endo"
 
pucker = "O4'-endo"
−
elif (phase < 144):
+
elif phase < 144:
 
pucker = "C1'-exo"
 
pucker = "C1'-exo"
−
elif (phase < 180):
+
elif phase < 180:
 
pucker = "C2'-endo"
 
pucker = "C2'-endo"
−
elif (phase < 216):
+
elif phase < 216:
 
pucker = "C3'-exo"
 
pucker = "C3'-exo"
−
elif (phase < 252):
+
elif phase < 252:
 
pucker = "C4'-endo"
 
pucker = "C4'-endo"
−
elif (phase < 288):
+
elif phase < 288:
 
pucker = "O4'-exo"
 
pucker = "O4'-exo"
−
elif (phase < 324):
+
elif phase < 324:
 
pucker = "C1'-endo"
 
pucker = "C1'-endo"
−
elif (phase < 360):
+
elif phase < 360:
 
pucker = "C2'-exo"
 
pucker = "C2'-exo"
 
else:
 
else:
Line 119: Line 194:
  
 
def store_atom(atom,residue):
 
def store_atom(atom,residue):
−
if (atom.name == "O4'" or atom.name == "O4*"):
+
if atom.name == "O4'" or atom.name == "O4*":
 
residue['O4*X'] = atom.coord[0]
 
residue['O4*X'] = atom.coord[0]
 
residue['O4*Y'] = atom.coord[1]
 
residue['O4*Y'] = atom.coord[1]
 
residue['O4*Z'] = atom.coord[2]
 
residue['O4*Z'] = atom.coord[2]
−
elif (atom.name == "C1'" or atom.name == "C1*"):
+
elif atom.name == "C1'" or atom.name == "C1*":
 
residue['C1*X'] = atom.coord[0]
 
residue['C1*X'] = atom.coord[0]
 
residue['C1*Y'] = atom.coord[1]
 
residue['C1*Y'] = atom.coord[1]
 
residue['C1*Z'] = atom.coord[2]
 
residue['C1*Z'] = atom.coord[2]
−
elif (atom.name == "C2'" or atom.name == "C2*"):
+
elif atom.name == "C2'" or atom.name == "C2*":
 
residue['C2*X'] = atom.coord[0]
 
residue['C2*X'] = atom.coord[0]
 
residue['C2*Y'] = atom.coord[1]
 
residue['C2*Y'] = atom.coord[1]
 
residue['C2*Z'] = atom.coord[2]
 
residue['C2*Z'] = atom.coord[2]
−
elif (atom.name == "C3'" or atom.name == "C3*"):
+
elif atom.name == "C3'" or atom.name == "C3*":
 
residue['C3*X'] = atom.coord[0]
 
residue['C3*X'] = atom.coord[0]
 
residue['C3*Y'] = atom.coord[1]
 
residue['C3*Y'] = atom.coord[1]
 
residue['C3*Z'] = atom.coord[2]
 
residue['C3*Z'] = atom.coord[2]
−
elif (atom.name == "C4'" or atom.name == "C4*"):
+
elif atom.name == "C4'" or atom.name == "C4*":
 
residue['C4*X'] = atom.coord[0]
 
residue['C4*X'] = atom.coord[0]
 
residue['C4*Y'] = atom.coord[1]
 
residue['C4*Y'] = atom.coord[1]
Line 150: Line 225:
 
for i in range (0,12):
 
for i in range (0,12):
 
C.append(residue[ribose[i]])
 
C.append(residue[ribose[i]])
−
theta['0']=dihedral(C)
+
theta['0']=calcdihedral(C)
 
 
 
C = []
 
C = []
Line 156: Line 231:
 
for i in range (0,12):
 
for i in range (0,12):
 
C.append(residue[ribose[i]])
 
C.append(residue[ribose[i]])
−
theta['1']=dihedral(C)
+
theta['1']=calcdihedral(C)
  
  
Line 163: Line 238:
 
for i in range (0,12):
 
for i in range (0,12):
 
C.append(residue[ribose[i]])
 
C.append(residue[ribose[i]])
−
theta['2']=dihedral(C)
+
theta['2']=calcdihedral(C)
  
  
Line 170: Line 245:
 
for i in range (0,12):
 
for i in range (0,12):
 
C.append(residue[ribose[i]])
 
C.append(residue[ribose[i]])
−
theta['3']=dihedral(C)
+
theta['3']=calcdihedral(C)
  
 
C = []
 
C = []
Line 176: Line 251:
 
for i in range (0,12):
 
for i in range (0,12):
 
C.append(residue[ribose[i]])
 
C.append(residue[ribose[i]])
−
theta['4']=dihedral(C)
+
theta['4']=calcdihedral(C)
 
 
 
return
 
return
Line 187: Line 262:
 
def shift_down(list,value):
 
def shift_down(list,value):
 
for i in range (0,value):
 
for i in range (0,value):
−
list.insert(len(list),list.pop(0))
+
list.append(list.pop(0))
 
return
 
return
  
−
def dihedral(C):
+
def calcdihedral(C):
 
 
 
DX12=C[0]-C[3]
 
DX12=C[0]-C[3]
Line 230: Line 305:
 
TS=1.0-DP12*DP12
 
TS=1.0-DP12*DP12
  
−
if (TS < 0):
+
if TS < 0:
 
TS=0
 
TS=0
 
else:
 
else:
Line 243: Line 318:
 
DP233=PX3*DX23+PY3*DY23+PZ3*DZ23
 
DP233=PX3*DX23+PY3*DY23+PZ3*DZ23
  
−
if (DP233 > 0):
+
if DP233 > 0:
 
ANGLE=-ANGLE
 
ANGLE=-ANGLE
  
 
ANGLE=math.degrees(ANGLE)
 
ANGLE=math.degrees(ANGLE)
  
−
if (ANGLE > 180):
+
if ANGLE > 180:
−
ANGLE=ANGLE-360
+
ANGLE-=360
−
if (ANGLE < -180):
+
if ANGLE < -180:
−
ANGLE=ANGLE+360
+
ANGLE+=360
  
 
return ANGLE
 
return ANGLE
  
 
cmd.extend("pucker",pucker)
 
cmd.extend("pucker",pucker)
 +
 
</source>
 
</source>
  
 
[[Category:Script_Library]]
 
[[Category:Script_Library]]
 +
[[Category:Biochemical_Properties]]
 +
[[Category:Biochemical_Scripts]]

Latest revision as of 03:36, 28 March 2014

Type Python Script
Download pucker.py
Author(s) Sean M. Law
License -
This code has been put under version control in the project Pymol-script-repo

DESCRIPTION

pucker.py is a PyMol script that returns the sugar pucker information (phase, amplitude, pucker) for a given selection.

This script uses its own dihedral calculation scheme rather than the get_dihedral command. Thus, it is lightning fast!

If a selection does not contain any ribose sugars then an error message is returned.

USAGE

load the script using the run command 

pucker selection
#Calculates the sugar information for the given selection

pucker selection, state=1
#Calculates the sugar information for the given selection in state 1

pucker selection, state=0 
#Iterates over all states and calculates the sugar information for the given selection

EXAMPLES

#fetch 1BNA.pdb
fetch 1bna

#Select DNA only
#Otherwise, you will get an error for water not having sugars
select DNA, not solvent

#Execute pucker command
pucker DNA

#The output should look like this
 Phase     Amp    Pucker  Residue
161.22   55.32  C2'-endo  A 1
139.52   41.67   C1'-exo  A 2
 92.82   38.31  O4'-endo  A 3
166.35   48.47  C2'-endo  A 4
128.57   46.30   C1'-exo  A 5
126.92   49.75   C1'-exo  A 6
101.30   47.32  O4'-endo  A 7
115.62   49.23   C1'-exo  A 8
140.44   46.37   C1'-exo  A 9
145.79   53.36  C2'-endo  A 10
147.47   47.04  C2'-endo  A 11
113.80   51.69   C1'-exo  A 12
 
 Phase     Amp    Pucker  Residue
153.24   43.15  C2'-endo  B 13
128.49   45.01   C1'-exo  B 14
 67.74   43.84   C4'-exo  B 15
149.33   41.01  C2'-endo  B 16
169.27   42.30  C2'-endo  B 17
147.03   42.30  C2'-endo  B 18
116.29   47.52   C1'-exo  B 19
129.62   49.92   C1'-exo  B 20
113.61   42.93   C1'-exo  B 21
156.34   50.98  C2'-endo  B 22
116.89   44.21   C1'-exo  B 23
 34.70   45.55  C3'-endo  B 24


PYMOL API

from pymol.cgo import *    # get constants
from math import *
from pymol import cmd

def pucker(selection,state=1):

	# Author: Sean Law
	# Institute: University of Michigan
	# E-mail: seanlaw@umich.edu

	#Comparison to output from 3DNA is identical
	#Note that the 3DNA output has the second chain
	#printed out in reverse order
	state=int(state)
	if state == 0:
		for state in range(1,cmd.count_states()+1):
			sel=cmd.get_model(selection,state)
			if state == 1:
				print " " #Blank line to separate chain output
				print "%5s  %6s  %6s  %8s  Residue" % ("State","Phase","Amp", "Pucker")
			get_pucker(sel,iterate=state)
	else:
		sel=cmd.get_model(selection,state)
		get_pucker(sel,iterate=0)
	return

def get_pucker (sel,iterate=0):
	iterate=int(iterate)
	first=1
	old=" "
	oldchain=" "
	residue={}
	theta={}
	count=0
	for atom in sel.atom:
		new=atom.chain+" "+str(atom.resi)
		newchain=atom.chain+" "+atom.segi
		if oldchain != newchain and first:
			if iterate == 0:
				print " " #Blank line to separate chain output
				print "%6s  %6s  %8s  Residue" % ("Phase", "Amp", "Pucker")
		if new != old and not first:
			#Check that all 5 atoms exist
			if len(residue) == 15:
				#Construct planes
				get_dihedrals(residue,theta)
				#Calculate pucker
				info = pseudo(residue,theta)
				print info+"  "+old
			else:
				print "There is no sugar in this residue"
			if oldchain != newchain:
				print " " #Blank line to separate chain output
				print "%6s  %6s  %8s  Residue" % ("Phase", "Amp", "Pucker")
			#Clear values
			residue={}
			dihedrals={}
			theta={}
			#Store new value
			store_atom(atom,residue)
		else:
			store_atom(atom,residue)
		first=0
		old=new
		oldchain=newchain
	#Final Residue
	#Calculate dihedrals for final residue
	if len(residue) == 15:
		#Construct planes
		get_dihedrals(residue,theta)
		#Calculate pucker for final residue
		info = pseudo(residue,theta)
		if iterate == 0:
			print info+"  "+old
		else:
			print "%5d  %s" % (iterate,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,theta):
	other=2*(sin(math.radians(36.0))+sin(math.radians(72.0)))
	
	#phase=atan2((theta4+theta1)-(theta3+theta0),2*theta2*(sin(math.radians(36.0))+sin(math.radians(72.0))))
	phase=atan2((theta['4']+theta['1'])-(theta['3']+theta['0']),theta['2']*other)
	amplitude=theta['2']/cos(phase)
	phase=math.degrees(phase)
	if phase < 0:
		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):
	if atom.name == "O4'" or atom.name == "O4*":
		residue['O4*X'] = atom.coord[0]
		residue['O4*Y'] = atom.coord[1]
		residue['O4*Z'] = atom.coord[2]
	elif atom.name == "C1'" or atom.name == "C1*":
		residue['C1*X'] = atom.coord[0]
		residue['C1*Y'] = atom.coord[1]
		residue['C1*Z'] = atom.coord[2]
	elif atom.name == "C2'" or atom.name == "C2*":
		residue['C2*X'] = atom.coord[0]
		residue['C2*Y'] = atom.coord[1]
		residue['C2*Z'] = atom.coord[2]
	elif atom.name == "C3'" or atom.name == "C3*":
		residue['C3*X'] = atom.coord[0]
		residue['C3*Y'] = atom.coord[1]
		residue['C3*Z'] = atom.coord[2]
	elif atom.name == "C4'" or atom.name == "C4*":
		residue['C4*X'] = atom.coord[0]
		residue['C4*Y'] = atom.coord[1]
		residue['C4*Z'] = atom.coord[2]
	return

def get_dihedrals(residue,theta):

	C = []
	ribose = residue.keys()
	ribose.sort()
	
	shift_up(ribose,6)
	for i in range (0,12):
		C.append(residue[ribose[i]])
	theta['0']=calcdihedral(C)
	
	C = []
	shift_down(ribose,3)
	for i in range (0,12):
		C.append(residue[ribose[i]])
	theta['1']=calcdihedral(C)


	C = []
	shift_down(ribose,3)
	for i in range (0,12):
		C.append(residue[ribose[i]])
	theta['2']=calcdihedral(C)


	C = []
	shift_down(ribose,3)
	for i in range (0,12):
		C.append(residue[ribose[i]])
	theta['3']=calcdihedral(C)

	C = []
	shift_down(ribose,3)
	for i in range (0,12):
		C.append(residue[ribose[i]])
	theta['4']=calcdihedral(C)
	
	return

def shift_up(list,value):
	for i in range (0,value):
		list.insert(0,list.pop())
	return

def shift_down(list,value):
	for i in range (0,value):
		list.append(list.pop(0))
	return

def calcdihedral(C):
	
	DX12=C[0]-C[3]
	DY12=C[1]-C[4]
	DZ12=C[2]-C[5]
	
	DX23=C[3]-C[6]
	DY23=C[4]-C[7]
	DZ23=C[5]-C[8]
	
	DX43=C[9]-C[6];
	DY43=C[10]-C[7];
	DZ43=C[11]-C[8];

	#Cross product to get normal
	
	PX1=DY12*DZ23-DY23*DZ12;
	PY1=DZ12*DX23-DZ23*DX12;
	PZ1=DX12*DY23-DX23*DY12;

	NP1=sqrt(PX1*PX1+PY1*PY1+PZ1*PZ1);

	PX1=PX1/NP1
	PY1=PY1/NP1
	PZ1=PZ1/NP1

	PX2=DY43*DZ23-DY23*DZ43;
	PY2=DZ43*DX23-DZ23*DX43;
	PZ2=DX43*DY23-DX23*DY43;

	NP2=sqrt(PX2*PX2+PY2*PY2+PZ2*PZ2);
	
	PX2=PX2/NP2
	PY2=PY2/NP2
	PZ2=PZ2/NP2

	DP12=PX1*PX2+PY1*PY2+PZ1*PZ2

	TS=1.0-DP12*DP12

	if TS < 0:
		TS=0
	else:
		TS=sqrt(TS)
	
	ANGLE=math.pi/2.0-atan2(DP12,TS)

	PX3=PY1*PZ2-PY2*PZ1
	PY3=PZ1*PX2-PZ2*PX1
	PZ3=PX1*PY2-PX2*PY1

	DP233=PX3*DX23+PY3*DY23+PZ3*DZ23

	if DP233 > 0:
		ANGLE=-ANGLE

	ANGLE=math.degrees(ANGLE)

	if ANGLE > 180:
		ANGLE-=360
	if ANGLE < -180:
		ANGLE+=360

	return ANGLE

cmd.extend("pucker",pucker)
Retrieved from "http://wiki.pymol.org/index.php?title=Pucker&oldid=9255"