AngleBetweenHelices - Revision history

Cchem: 1 revision

2017-01-09T05:54:53Z

1 revision

Speleo3: psico

2016-12-15T14:38:45Z

psico

Pyadmin: 13 revisions

2014-03-28T01:06:19Z

13 revisions

Tlinnet at 15:00, 15 January 2012

2012-01-15T15:00:50Z

Speleo3: gallery float right

2011-12-11T09:53:37Z

gallery float right

Speleo3: Infobox script-repo

2011-12-09T09:41:10Z

Infobox script-repo

Speleo3: use Template:PymolScriptRepoDownload

2011-12-04T11:53:40Z

use Template:PymolScriptRepoDownload

Tlinnet: /* The Script */

2011-12-03T17:48:10Z

The Script

Show changes

Tlinnet: /* Example */

2011-12-03T17:46:09Z

Example

New page

<gallery perrow=2 heights=191px widths=250px>
Image:Helix_orientation1.png|Orientation of a helix
Image:Helix_orientation2.png|Orientation of two helices, their angle that separates the two is 145.08 degrees as reported by the script
</gallery>
Calculate angle between alpha-helices or beta-sheets.
There are four different methods implemented to fit a helix, two of them also work for sheets or loops.

= Commands =
<source lang="python">
helix_orientation selection [, visualize [, sigma_cutoff [, quiet ]]]

helix_orientation_hbond selection [, visualize [, cutoff [, quiet ]]]

loop_orientation selection [, visualize [, quiet ]]

cafit_orientation selection [, visualize [, quiet ]]

angle_between_helices selection1, selection2 [, method [, visualize [, quiet ]]]
</source>

= Example =

<syntaxhighlight lang="python">
import anglebetweenhelices

fetch 2x19, async=0
select hel1, /2x19//B/23-36/
select hel2, /2x19//B/40-54/

# just calculate/visualize orientation of single alpha-helix
helix_orientation_hbond hel1

# get angle between two helices
angle_between_helices hel1, hel2
angle_between_helices hel1, hel2, method=1
angle_between_helices hel1, hel2, method=2

# get angle between beta-sheets
select sheet1, A/47-54/
select sheet6, A/146-149/
angle_between_helices sheet1, sheet6, method=loop_orientation
angle_between_helices sheet1, sheet6, method=cafit_orientation
</syntaxhighlight>

Output:
<syntaxhighlight lang="python">
PyMOL>angle_between_helices hel1, hel2, method=cafit_orientation
Using method: cafit_orientation
Angle: 145.08 deg
</syntaxhighlight>

= The Script =
<source lang="python">
'''
(c) 2010 Thomas Holder
'''

from pymol import cmd, stored, CmdException
from chempy import cpv
import math

def _vec_sum(vec_list):
# this is the same as
# return numpy.array(vec_list).sum(0).tolist()
vec = cpv.get_null()
for x in vec_list:
vec = cpv.add(vec, x)
return vec

def _mean_and_std(x):
# this is the same as
# return (numpy.mean(x), numpy.std(x, ddof=1))
N = len(x)
if N < 2:
return (x[0], 0.0)
mu = sum(x)/float(N)
var = sum((i - mu)**2 for i in x) / float(N - 1)
return (mu, var**0.5)

def _common_orientation(selection, vec, visualize=1, quiet=0):
'''
Common part of different helix orientation functions. Does calculate
the center of mass and does the visual feedback.
'''
stored.x = []
cmd.iterate_state(-1, '(%s) and name CA' % (selection),
'stored.x.append([x,y,z])')
if len(stored.x) < 2:
print 'warning: count(CA) < 2'
raise CmdException
center = cpv.scale(_vec_sum(stored.x), 1./len(stored.x))
if visualize:
scale = cpv.distance(stored.x[0], stored.x[-1])
visualize_orientation(vec, center, scale, True)
cmd.zoom(selection, buffer=2)
if not quiet:
print 'Center: (%.2f, %.2f, %.2f) Direction: (%.2f, %.2f, %.2f)' % tuple(center + vec)
return center, vec

def visualize_orientation(direction, center=[0,0,0], scale=1.0, symmetric=False, color='green', color2='red'):
'''
Draw an arrow. Helper function for "helix_orientation" etc.
'''
from pymol import cgo
color_list = cmd.get_color_tuple(color)
color2_list = cmd.get_color_tuple(color2)
if symmetric:
scale *= 0.5
end = cpv.add(center, cpv.scale(direction, scale))
radius = 0.3
obj = [cgo.SAUSAGE]
obj.extend(center)
obj.extend(end)
obj.extend([
radius,
0.8, 0.8, 0.8,
])
obj.extend(color_list)
if symmetric:
start = cpv.sub(center, cpv.scale(direction, scale))
obj.append(cgo.SAUSAGE)
obj.extend(center)
obj.extend(start)
obj.extend([
radius,
0.8, 0.8, 0.8,
])
obj.extend(color2_list)
coneend = cpv.add(end, cpv.scale(direction, 4.0*radius))
obj.append(cgo.CONE)
obj.extend(end)
obj.extend(coneend)
obj.extend([
radius * 1.75,
0.0,
])
obj.extend(color_list * 2)
obj.extend([
1.0, 1.0, # Caps
])
cmd.load_cgo(obj, cmd.get_unused_name('oriVec'), zoom=0)

def cafit_orientation(selection, visualize=1, quiet=0):
'''
DESCRIPTION

Get the center and direction of a peptide by least squares
linear fit on CA atoms.

USAGE

cafit_orientation selection [, visualize]

NOTES

Requires python module "numpy".

SEE ALSO

helix_orientation
'''
visualize, quiet = int(visualize), int(quiet)
import numpy
stored.x = list()
cmd.iterate_state(-1, '(%s) and name CA' % (selection),
'stored.x.append([x,y,z])')
x = numpy.array(stored.x)
U,s,Vh = numpy.linalg.svd(x - x.mean(0))
vec = cpv.normalize(Vh[0])
if cpv.dot_product(vec, x[-1] - x[0]) < 0:
vec = cpv.negate(vec)
return _common_orientation(selection, vec, visualize, quiet)

def loop_orientation(selection, visualize=1, quiet=0):
'''
DESCRIPTION

Get the center and approximate direction of a peptide. Works for any
secondary structure.
Averages direction of N(i)->C(i) pseudo bonds.

USAGE

loop_orientation selection [, visualize]

SEE ALSO

helix_orientation
'''
visualize, quiet = int(visualize), int(quiet)
stored.x = dict()
cmd.iterate_state(-1, '(%s) and name N+C' % (selection),
'stored.x.setdefault(chain + resi, dict())[name] = x,y,z')
vec = cpv.get_null()
count = 0
for x in stored.x.itervalues():
if 'C' in x and 'N' in x:
vec = cpv.add(vec, cpv.sub(x['C'], x['N']))
count += 1
if count == 0:
print 'warning: count == 0'
raise CmdException
vec = cpv.normalize(vec)
return _common_orientation(selection, vec, visualize, quiet)

def helix_orientation(selection, visualize=1, sigma_cutoff=1.5, quiet=0):
'''
DESCRIPTION

Get the center and direction of a helix as vectors. Will only work
for helices and gives slightly different results than loop_orientation.
Averages direction of C(i)->O(i) bonds.

USAGE

helix_orientation selection [, visualize [, sigma_cutoff]]

ARGUMENTS

selection = string: atom selection of helix

visualize = 0 or 1: show fitted vector as arrow {default: 1}

sigma_cutoff = float: drop outliers outside
(standard_deviation * sigma_cutoff) {default: 1.5}

SEE ALSO

angle_between_helices, helix_orientation_hbond, loop_orientation, cafit_orientation
'''
visualize, quiet, sigma_cutoff = int(visualize), int(quiet), float(sigma_cutoff)
stored.x = dict()
cmd.iterate_state(-1, '(%s) and name C+O' % (selection),
'stored.x.setdefault(chain + resi, dict())[name] = x,y,z')
vec_list = []
count = 0
for x in stored.x.itervalues():
if 'C' in x and 'O' in x:
vec_list.append(cpv.sub(x['O'], x['C']))
count += 1
if count == 0:
print 'warning: count == 0'
raise CmdException
vec = _vec_sum(vec_list)
if count > 2 and sigma_cutoff > 0:
angle_list = [cpv.get_angle(vec, x) for x in vec_list]
angle_mu, angle_sigma = _mean_and_std(angle_list)
vec_list = [vec_list[i] for i in range(len(vec_list))
if abs(angle_list[i] - angle_mu) < angle_sigma * sigma_cutoff]
if not quiet:
print 'Dropping %d outlier(s)' % (len(angle_list) - len(vec_list))
vec = _vec_sum(vec_list)
vec = cpv.normalize(vec)
return _common_orientation(selection, vec, visualize, quiet)

def helix_orientation_hbond(selection, visualize=1, cutoff=3.5, quiet=0):
'''
DESCRIPTION

Get the center and direction of a helix as vectors. Will only work
for alpha helices and gives slightly different results than
helix_orientation. Averages direction of O(i)->N(i+4) hydrogen bonds.

USAGE

helix_orientation selection [, visualize [, cutoff]]

ARGUMENTS

cutoff = float: maximal hydrogen bond distance {default: 3.5}

SEE ALSO

helix_orientation
'''
visualize, quiet, cutoff = int(visualize), int(quiet), float(cutoff)
stored.x = dict()
cmd.iterate_state(-1, '(%s) and name N+O' % (selection),
'stored.x.setdefault(resv, dict())[name] = x,y,z')
vec_list = []
for resi in stored.x:
resi_other = resi + 4
if 'O' in stored.x[resi] and resi_other in stored.x:
if 'N' in stored.x[resi_other]:
vec = cpv.sub(stored.x[resi_other]['N'], stored.x[resi]['O'])
if cpv.length(vec) < cutoff:
vec_list.append(vec)
if len(vec_list) == 0:
print 'warning: count == 0'
raise CmdException
vec = _vec_sum(vec_list)
vec = cpv.normalize(vec)
return _common_orientation(selection, vec, visualize, quiet)

def angle_between_helices(selection1, selection2, method='helix_orientation', visualize=1, quiet=0):
'''
DESCRIPTION

Calculates the angle between two helices

USAGE

angle_between_helices selection1, selection2 [, method [, visualize]]

ARGUMENTS

selection1 = string: atom selection of first helix

selection2 = string: atom selection of second helix

method = string: function to calculate orientation {default: helix_orientation}
or int: 0: helix_orientation, 1: helix_orientation_hbond,
2: loop_orientation, 3: cafit_orientation

visualize = 0 or 1: show fitted vector as arrow {default: 1}

SEE ALSO

helix_orientation, helix_orientation_hbond, loop_orientation, cafit_orientation
'''
visualize, quiet = int(visualize), int(quiet)
methods = {
'0': helix_orientation,
'1': helix_orientation_hbond,
'2': loop_orientation,
'3': cafit_orientation,
}
methods.update((x.__name__, x) for x in methods.values())
try:
orientation = methods[str(method)]
except KeyError:
print 'no such method:', method
raise CmdException
if not quiet:
print 'Using method:', orientation.__name__
cen1, dir1 = orientation(selection1, visualize, quiet=1)
cen2, dir2 = orientation(selection2, visualize, quiet=1)
angle = cpv.get_angle(dir1, dir2)
angle_deg = math.degrees(angle)
if not quiet:
print 'Angle: %.2f deg' % (angle_deg)
if visualize:
cmd.zoom('(%s) or (%s)' % (selection1, selection2), buffer=2)
return angle_deg

cmd.extend('helix_orientation', helix_orientation)
cmd.extend('helix_orientation_hbond', helix_orientation_hbond)
cmd.extend('loop_orientation', loop_orientation)
cmd.extend('cafit_orientation', cafit_orientation)
cmd.extend('angle_between_helices', angle_between_helices)
</source>

[[Category:Plugins]]
[[Category:Structural_Biology_Scripts]]

Speleo3: new page

2010-09-13T09:16:35Z

new page

New page

Calculate angle between alpha-helices or beta-sheets.
There are four different methods implemented to fit a helix, two of them also work for sheets or loops.

= Commands =

* helix_orientation selection [, visualize [, sigma_cutoff [, quiet ]]]
* helix_orientation_hbond selection [, visualize [, cutoff [, quiet ]]]
* loop_orientation selection [, visualize [, quiet ]]
* cafit_orientation selection [, visualize [, quiet ]]
* angle_between_helices selection1, selection2 [, method [, visualize [, quiet ]]]

= Example =

<source lang="python">
fetch 2x19
select hel1, /2x19//B/23-36/
select hel2, /2x19//B/40-54/

# just calculate/visualize orientation of single alpha-helix
helix_orientation_hbond hel1

# get angle between two helices
angle_between_helices hel1, hel2
angle_between_helices hel1, hel2, method=1
angle_between_helices hel1, hel2, method=2

# get angle between beta-sheets
select sheet1, A/47-54/
select sheet6, A/146-149/
angle_between_helices sheet1, sheet6, method=loop_orientation
angle_between_helices sheet1, sheet6, method=cafit_orientation
</source>

= The Script =
<source lang="python">
'''
(c) 2010 Thomas Holder
'''

from pymol import cmd, stored, CmdException
from chempy import cpv
import math

def _vec_sum(vec_list):
# this is the same as
# return numpy.array(vec_list).sum(0).tolist()
vec = cpv.get_null()
for x in vec_list:
vec = cpv.add(vec, x)
return vec

def _mean_and_std(x):
# this is the same as
# return (numpy.mean(x), numpy.std(x, ddof=1))
N = len(x)
if N < 2:
return (x[0], 0.0)
mu = sum(x)/float(N)
var = sum((i - mu)**2 for i in x) / float(N - 1)
return (mu, var**0.5)

def _common_orientation(selection, vec, visualize=1, quiet=0):
'''
Common part of different helix orientation functions. Does calculate
the center of mass and does the visual feedback.
'''
stored.x = []
cmd.iterate_state(-1, '(%s) and name CA' % (selection),
'stored.x.append([x,y,z])')
if len(stored.x) < 2:
print 'warning: count(CA) < 2'
raise CmdException
center = cpv.scale(_vec_sum(stored.x), 1./len(stored.x))
if visualize:
scale = cpv.distance(stored.x[0], stored.x[-1])
visualize_orientation(vec, center, scale, True)
cmd.zoom(selection, buffer=2)
if not quiet:
print 'Center: (%.2f, %.2f, %.2f) Direction: (%.2f, %.2f, %.2f)' % tuple(center + vec)
return center, vec

def visualize_orientation(direction, center=[0,0,0], scale=1.0, symmetric=False, color='green', color2='red'):
'''
Draw an arrow. Helper function for "helix_orientation" etc.
'''
from pymol import cgo
color_list = cmd.get_color_tuple(color)
color2_list = cmd.get_color_tuple(color2)
if symmetric:
scale *= 0.5
end = cpv.add(center, cpv.scale(direction, scale))
radius = 0.3
obj = [cgo.SAUSAGE]
obj.extend(center)
obj.extend(end)
obj.extend([
radius,
0.8, 0.8, 0.8,
])
obj.extend(color_list)
if symmetric:
start = cpv.sub(center, cpv.scale(direction, scale))
obj.append(cgo.SAUSAGE)
obj.extend(center)
obj.extend(start)
obj.extend([
radius,
0.8, 0.8, 0.8,
])
obj.extend(color2_list)
coneend = cpv.add(end, cpv.scale(direction, 4.0*radius))
obj.append(cgo.CONE)
obj.extend(end)
obj.extend(coneend)
obj.extend([
radius * 1.75,
0.0,
])
obj.extend(color_list * 2)
obj.extend([
1.0, 1.0, # Caps
])
cmd.load_cgo(obj, cmd.get_unused_name('oriVec'), zoom=0)

def cafit_orientation(selection, visualize=1, quiet=0):
'''
DESCRIPTION

Get the center and direction of a peptide by least squares
linear fit on CA atoms.

USAGE

cafit_orientation selection [, visualize]

NOTES

Requires python module "numpy".

SEE ALSO

helix_orientation
'''
visualize, quiet = int(visualize), int(quiet)
import numpy
stored.x = list()
cmd.iterate_state(-1, '(%s) and name CA' % (selection),
'stored.x.append([x,y,z])')
x = numpy.array(stored.x)
U,s,Vh = numpy.linalg.svd(x - x.mean(0))
vec = cpv.normalize(Vh[0])
if cpv.dot_product(vec, x[-1] - x[0]) < 0:
vec = cpv.negate(vec)
return _common_orientation(selection, vec, visualize, quiet)

def loop_orientation(selection, visualize=1, quiet=0):
'''
DESCRIPTION

Get the center and approximate direction of a peptide. Works for any
secondary structure.
Averages direction of N(i)->C(i) pseudo bonds.

USAGE

loop_orientation selection [, visualize]

SEE ALSO

helix_orientation
'''
visualize, quiet = int(visualize), int(quiet)
stored.x = dict()
cmd.iterate_state(-1, '(%s) and name N+C' % (selection),
'stored.x.setdefault(chain + resi, dict())[name] = x,y,z')
vec = cpv.get_null()
count = 0
for x in stored.x.itervalues():
if 'C' in x and 'N' in x:
vec = cpv.add(vec, cpv.sub(x['C'], x['N']))
count += 1
if count == 0:
print 'warning: count == 0'
raise CmdException
vec = cpv.normalize(vec)
return _common_orientation(selection, vec, visualize, quiet)

def helix_orientation(selection, visualize=1, sigma_cutoff=1.5, quiet=0):
'''
DESCRIPTION

Get the center and direction of a helix as vectors. Will only work
for helices and gives slightly different results than loop_orientation.
Averages direction of C(i)->O(i) bonds.

USAGE

helix_orientation selection [, visualize [, sigma_cutoff]]

ARGUMENTS

selection = string: atom selection of helix

visualize = 0 or 1: show fitted vector as arrow {default: 1}

sigma_cutoff = float: drop outliers outside
(standard_deviation * sigma_cutoff) {default: 1.5}

SEE ALSO

angle_between_helices, helix_orientation_hbond, loop_orientation, cafit_orientation
'''
visualize, quiet, sigma_cutoff = int(visualize), int(quiet), float(sigma_cutoff)
stored.x = dict()
cmd.iterate_state(-1, '(%s) and name C+O' % (selection),
'stored.x.setdefault(chain + resi, dict())[name] = x,y,z')
vec_list = []
count = 0
for x in stored.x.itervalues():
if 'C' in x and 'O' in x:
vec_list.append(cpv.sub(x['O'], x['C']))
count += 1
if count == 0:
print 'warning: count == 0'
raise CmdException
vec = _vec_sum(vec_list)
if count > 2 and sigma_cutoff > 0:
angle_list = [cpv.get_angle(vec, x) for x in vec_list]
angle_mu, angle_sigma = _mean_and_std(angle_list)
vec_list = [vec_list[i] for i in range(len(vec_list))
if abs(angle_list[i] - angle_mu) < angle_sigma * sigma_cutoff]
if not quiet:
print 'Dropping %d outlier(s)' % (len(angle_list) - len(vec_list))
vec = _vec_sum(vec_list)
vec = cpv.normalize(vec)
return _common_orientation(selection, vec, visualize, quiet)

def helix_orientation_hbond(selection, visualize=1, cutoff=3.5, quiet=0):
'''
DESCRIPTION

Get the center and direction of a helix as vectors. Will only work
for alpha helices and gives slightly different results than
helix_orientation. Averages direction of O(i)->N(i+4) hydrogen bonds.

USAGE

helix_orientation selection [, visualize [, cutoff]]

ARGUMENTS

cutoff = float: maximal hydrogen bond distance {default: 3.5}

SEE ALSO

helix_orientation
'''
visualize, quiet, cutoff = int(visualize), int(quiet), float(cutoff)
stored.x = dict()
cmd.iterate_state(-1, '(%s) and name N+O' % (selection),
'stored.x.setdefault(resv, dict())[name] = x,y,z')
vec_list = []
for resi in stored.x:
resi_other = resi + 4
if 'O' in stored.x[resi] and resi_other in stored.x:
if 'N' in stored.x[resi_other]:
vec = cpv.sub(stored.x[resi_other]['N'], stored.x[resi]['O'])
if cpv.length(vec) < cutoff:
vec_list.append(vec)
if len(vec_list) == 0:
print 'warning: count == 0'
raise CmdException
vec = _vec_sum(vec_list)
vec = cpv.normalize(vec)
return _common_orientation(selection, vec, visualize, quiet)

def angle_between_helices(selection1, selection2, method='helix_orientation', visualize=1, quiet=0):
'''
DESCRIPTION

Calculates the angle between two helices

USAGE

angle_between_helices selection1, selection2 [, method [, visualize]]

ARGUMENTS

selection1 = string: atom selection of first helix

selection2 = string: atom selection of second helix

method = string: function to calculate orientation {default: helix_orientation}
or int: 0: helix_orientation, 1: helix_orientation_hbond,
2: loop_orientation, 3: cafit_orientation

visualize = 0 or 1: show fitted vector as arrow {default: 1}

SEE ALSO

helix_orientation, helix_orientation_hbond, loop_orientation, cafit_orientation
'''
visualize, quiet = int(visualize), int(quiet)
methods = {
'0': helix_orientation,
'1': helix_orientation_hbond,
'2': loop_orientation,
'3': cafit_orientation,
}
methods.update((x.__name__, x) for x in methods.values())
try:
orientation = methods[str(method)]
except KeyError:
print 'no such method:', method
raise CmdException
if not quiet:
print 'Using method:', orientation.__name__
cen1, dir1 = orientation(selection1, visualize, quiet=1)
cen2, dir2 = orientation(selection2, visualize, quiet=1)
angle = cpv.get_angle(dir1, dir2)
angle_deg = math.degrees(angle)
if not quiet:
print 'Angle: %.2f deg' % (angle_deg)
if visualize:
cmd.zoom('(%s) or (%s)' % (selection1, selection2), buffer=2)
return angle_deg

cmd.extend('helix_orientation', helix_orientation)
cmd.extend('helix_orientation_hbond', helix_orientation_hbond)
cmd.extend('loop_orientation', loop_orientation)
cmd.extend('cafit_orientation', cafit_orientation)
cmd.extend('angle_between_helices', angle_between_helices)
</source>

[[Category:Plugins]]
[[Category:Structural_Biology_Scripts]]

@@ Line 4: / Line 4: @@
 |author    = [[User:Speleo3|Thomas Holder]]
 |license   = BSD
 }}
@@ Line 15: / Line 19: @@
 = Commands =
-helix_orientation_hbond selection [, visualize [, cutoff [, quiet ]]]
+ angle_between_helices selection1, selection2 [, method [, visualize ]]
-loop_orientation selection [, visualize [, quiet ]]
+ helix_orientation_hbond selection [, visualize [, cutoff ]]
-cafit_orientation selection [, visualize [, quiet ]]
+ loop_orientation selection [, visualize ]
-angle_between_helices selection1, selection2 [, method [, visualize [, quiet ]]]
+ cafit_orientation selection [, visualize ]
 = Example =
@@ Line 57: / Line 60: @@
 Angle: 145.08 deg
 </syntaxhighlight>
 [[Category:Structural_Biology_Scripts]]
 [[Category:Pymol-script-repo]]

← Older revision		Revision as of 15:00, 15 January 2012
Line 59:		Line 59:

	[[Category:Structural_Biology_Scripts]]		[[Category:Structural_Biology_Scripts]]
		+	[[Category:Pymol-script-repo]]

← Older revision	Revision as of 05:54, 9 January 2017
(No difference)

← Older revision	Revision as of 01:06, 28 March 2014
(No difference)

@@ Line 6: / Line 6: @@
 }}
-<gallery perrow=2 heights=191px widths=250px>
+<gallery perrow=1 heights=191px widths=250px style="float:right; clear: right">
 Image:Helix_orientation1.png|Orientation of a helix
 Image:Helix_orientation2.png|Orientation of two helices, their angle that separates the two is 145.08 degrees as reported by the script
 </gallery>
 Calculate angle between alpha-helices or beta-sheets.
 There are four different methods implemented to fit a helix, two of them also work for sheets or loops.

@@ Line 51: / Line 51: @@
 == Python Code ==
 [[Category:Plugins]]
 [[Category:Structural_Biology_Scripts]]