Python Programming for Biology: Bioinformatics and Beyond

Figure 9.3.  A cartoon visualisation of the torsion angle defined by a sequence of

four points. A torsion angle, also known as a dihedral angle, is the measurement of the

amount of twist between two planes about the axis of their intersection. In three

dimensions a torsion angle may be defined by four distinct points, as is commonly found

in molecular structures, where the points represent atoms, and where the central two atoms

define the intersection axis and each end atom defines the direction of a plane.

Inside  the  function  the  first  procedure  is  to  define  three  vectors,  which  may  represent

the  bonds  between  three  pairs  of  atom  coordinates,  by  subtracting  coordinates  from  one

another; the resulting difference vector represents how to move from one coordinate to the

other. Two perpendicular cross-product vectors are then calculated: for the first and second

---

vectors and for the second and third vectors. The basic idea here is that even though the

cross-product  vectors  are  perpendicular  to  the  second  (central)  vector  they  will  have  the

same  angle  of  twist  (relative  to  the  central  vector)  between  them  as  the  first  and  third

vectors  do.  However,  because  they  are  perpendicular  to  the  central  vector  the  angle

between them is the angle we want: the angle of twist around the axis of the central vector.

Hence,  the  perpendicular  vectors  are  then  used  to  calculate  three  dot  products.  The  first

represents the ‘shadow’ projection of one vector on the other, which naturally depends on

the angle between them, and the other two dot products will simply be the lengths of the

vectors squared. From these the cosine of the torsion angle is obtained.

def calcTorsionAngle(coord1, coord2, coord3, coord4):

bondVec12 = coord1 - coord2

bondVec32 = coord3 - coord2

bondVec43 = coord4 - coord3

perpVec13 = cross(bondVec12, bondVec32)

perpVec24 = cross(bondVec43, bondVec32)

projection = dot(perpVec13, perpVec24)

squareDist13 = dot(perpVec13, perpVec13)

squareDist24 = dot(perpVec24, perpVec24)

cosine = projection / sqrt(squareDist13*squareDist24)

cosine = min(1.0, max(-1.0, cosine))

angle = acos(cosine)

if dot(perpVec13, cross(perpVec24, bondVec32)) < 0:

angle = -angle

return angle

Note  that  we  check  to  make  sure  the  cosine  is  between  −1  and  +1;  we  might  have  a

small error in the calculation from the use of floating point numbers so we use the min()

and  max()  functions  to  reinforce  the  limits.  The  angle  is  calculated  by  using  the  inverse

cosine  function  acos()  (meaning  arccosine).  Finally,  we  use  a  dot  and  cross-product  to

check  whether  we  actually  calculated  a  negative  angle;  because  the  inverse  cosine  will

only give a value between zero and π radians (180°) we don’t otherwise know whether the

first coordinate is ‘above’ or ‘below’ the plane defined by the other three. We can then test

the function:

from numpy import array, degrees

p1 = array([2,1,1])

p2 = array([2,0,0])

p3 = array([3,0,0])

p4 = array([3,1,-1])

angle = calcTorsionAngle(p1, p2, p3, p4)

print( degrees(angle) ) # -90.0

1

 NumPy  relies  on  fast  calculation  routines  from  a  library  called  LAPACK

(

http://www.netlib.org/lapack/

).

---