Python Programming for Biology: Bioinformatics and Beyond

Distances and angles

The  next  examples  involve  looking  inside  structures  to  extract  particular  measurements,

like atomic distances and angles. Specifically, the angles we will be looking at are torsion

(or dihedral) angles, which represent the angle of twist about an axis; here the axis will be

along a chemical bond and the twist will be the relative orientation of two neighbouring

bonds.

Firstly,  to  get  the  distances  between  two  atoms  is  very  easy:  we  simply  get  the

coordinates  of  each  atom,  get  the  difference  between  the  coordinates  (stored  in  the

variable deltas), add together the squares of these differences, and finally take the square

root  of  the  total.

12

 Naturally  because  we  have  NumPy  arrays  most  of  this  is  done  in  an

element-wise manner. It should be noted that the sum of squares can be calculated as the

dot  product  of  deltas  with  itself  (though  we  could  also  do  (deltas*deltas).sum(axis=1)).

Because  we  are  using  dot  on  one-dimensional  arrays  (vectors)  in  this  instance  we  get  a

single number out, i.e. we are not using it for matrix multiplication. So given two of our

Structure.Atom objects, named atomA and atomB, we do:

from math import sqrt

deltas = atomA.coords – atomB.coords

sumSquares = dot(deltas, deltas)

distance = sqrt( sumSquares )

---

We move on from this toy example to a function that takes a structure and gives back a

list of atoms that are within a certain distance from a specified point. This might be useful

if you wanted to work out which other atoms a given atom may be interacting with. The

function  takes  a  structure  object,  position  and  a  distance  limit,  to  specify  the  catchment

radius relative to the position. Initially within the function a blank list is created to collect

the close atoms, then we convert the input xyz location into a NumPy array and we extract

the  coordinate  array  from  the  structure  object  as  illustrated  previously.  Next  we  define

limit2 as the square of the limiting search radius. This is done so that when we are looking

though the coordinates we do not have to repeatedly calculate square roots, which would

needlessly  slow  things  down.  Given  that  distances  are  always  positive  numbers,  we  can

say whether one distance is larger than another (the search limit in this case) by comparing

the squares of the distances, rather than the actual distances.

def findCloseAtoms(structure, xyz, limit=5.0):

closeAtoms = []

xyz = array(xyz)

atoms, coords = getAtomCoords(structure)

limit2 = limit * limit

deltas = coords – xyz

squares = deltas * deltas

sumSquares = squares.sum(axis=1)

boolArray = sumSquares < limit2

indices = boolArray.nonzero()[0] # Array for first and only axis

closeAtoms = [atoms[i] for i in indices]

return closeAtoms

Next in the function deltas is calculated as the difference from the atom coordinates to

the  xyz  position.  The  square  distance  for  each  atom’s  separation  is  calculated  by

multiplying  the  differences  with  themselves  and  then  summing  along  the  spatial

dimensions  (axis=1)  to  give  sumSquares.  These  are  then  compared  to  the  square  of  the

limiting radius. Using comparisons on NumPy arrays will act in an element-wise manner

to  generate  an  array  of  Boolean  values.  Hence  boolArray  will  have  true  elements  where

the square distance is less than the limit. The indices of the non-zero (true) elements are

extracted with .nonzero(), taking note that this gives a tuple with separate indices for each

array  axis.  Thus,  because  boolArray  is  only  one-dimensional,  here  we  take  the  first  and

only  array  of  indices  ([0]).  These  indices  can  finally  be  used  to  select  the  appropriate

group of close atoms from the larger list via a list comprehension. This list of close atoms

is  then  passed  back  at  the  end  of  the  function.  We  can  test  the  function  with  a  Structure

object,  here  finding  atoms  that  are  close  to  the  centre  of  mass  coordinates  that  were

determined earlier:

atoms = findCloseAtoms(struc, (18.89, 0.36, 1.24), 5.0)

for atom in atoms:

print(atom.residue.code, atom.name)

Now we will move on to calculate the twist angle between atoms, imagined as a group

---

of  four  atoms  in  a  rough  rectangle,  connected  by  three  bonds  so  that  three  edges  are

connected  and  one  edge  is  not.  The  angle  we  will  be  measuring  is  the  amount  of  twist

between  the  first  and  last  bonds,  around  the  central  bond.  This  can  be  visualised  by

looking directly along the central bond, so that the other two bonds make a ‘V’ shape; the

torsion angle is the angle made by this shape.

We can use the function described in

Chapter 9

that calculates torsion angles to measure

the

ϕ and ψ angles for an amino acid residue we find in a protein structure. These angles

represent  the  twist  of  the  polypeptide  chain  backbone  either  side  of  the  residue’s  alpha-

carbon  atom  position  (where  the  side  chain  comes  out).  The  input  to  the  function  is

naturally  a  Residue  object,  from  our  structure  data  model,  and  an  option  to  specify

whether we want the angle to be represented in units of degrees (the default) or in radians.

Because of the choice of angle units, we also import a standard function which allows us

to do the conversion between degrees and radians.

13

from Maths import calcTorsionAngle

from math import degrees

def getPhiPsi(residue, inDegrees=True):

Initially  the

ϕ  and  ψ  angles  are  initialised  as  None,  so  that  if  we  fail  to  find  certain

backbone  atoms  (this  happens  for  the  first  and  last  residue  of  a  chain)  there  is  an

indication that the angles were not definable. Next we get a list of all of the residues in the

input residue’s chain, so that we can get hold of the residues either side. Given that we are

using our data model objects, getting the list of residues is a simple matter of first finding

the  Chain  object  to  which  the  input  Residue  object  belongs  (this  is  the  parent  link),  and

then following the .residues link to get the list of all residues (children) for the chain.

phi = None

psi = None

chain = residue.chain

residues = chain.residues

For  the  input  residue  we  find  three  Atom  objects  representing  the  polypeptide

backbone,  by  fetching  atoms  with  the  appropriate  names.  For  each  atom  we  collect  the

array of their coordinates.

atomN = residue.getAtom('N')

atomCa = residue.getAtom('CA')

atomC = residue.getAtom('C')

coordsN = atomN.coords

coordsCa = atomCa.coords

coordsC = atomC.coords

A  positional  index  is  defined  to  represent  the  location  of  the  current  residue  in  the

polypeptide chain.

index = residues.index(residue)

---

if index > 0:

An index greater than zero means we are not at the start and thus can find the previous

residue in the chain, by taking one from the index and fetching the residue from the list.

We use the previous residue to fetch the atom, and then the coordinates, for the preceding

carbonyl carbon (‘C’) atom, which is the last of the positions required to calculate the phi

angle. Armed with the four coordinates, and using them in the correct order, the angle is

calculated using the function described previously. Then the angle is converted to units of

degrees if required.

residuePrev = residues[index-1]

atomC0 = residuePrev.getAtom('C')

coordsC0 = atomC0.coords

phi = calcTorsionAngle(coordsC0, coordsN, coordsCa, coordsC)

if inDegrees:

phi = degrees(phi)

If  the  index  is  less  than  the  last  possible  position,  then  we  are  able  to  find  the  next

residue in the chain, by adding one to the index, and thus also are able to calculate the psi

angle  in  the  same  manner.  This  torsion  angle  requires  the  coordinates  of  the  backbone

nitrogen  atom  in  the  next  residue,  i.e.  it  is  on  the  opposite  side  of  the  input  residue,

compared to phi.

if index < ( len(residues)-1 ):

residueNext = residues[index+1]

atomN2 = residueNext.getAtom('N')

coordsN2 = atomN2.coords

psi = calcTorsionAngle(coordsN, coordsCa, coordsC, coordsN2)

if inDegrees:

psi = degrees(psi)

return phi, psi

Finally, the pair of angles is returned from the function. We can test the code by calling

the  function  on  all  of  the  residues  in  a  structure  (although  strictly  speaking  we  ought  to

check that it is a protein first). Note that we exclude the first and last residues by using the

slice notation [1:-1]. The measured

ϕ and ψ angles are placed into separate lists, which are

then passed as input to a function from the pyplot module to generate a scatter plot which

allows  us  to  visualise  the  data  (see

http://www.cambridge.org/pythonforbiology

 for

Matplotlib download sites). This kind of display is called a Ramachandran plot, after one

of  its  inventors,  and  will  often  have  background  colours  indicating  the  likelihood  of

finding  different  angles.  Positions  of  the

ϕ  and  ψ  angles  on  the  Ramachandran  plot  are

indicative  of  the  secondary  structure  of  the  residues.  Although  each  secondary-structure

type (defined in terms of hydrogen bonding) will cover a range of angle values, a typical

alpha helix will be at (

ϕ = −60°, ψ = −45°) and a typical beta-strand at (ϕ = −135°, ψ =

---

135°).

from matplotlib import pyplot

phiList = []

psiList = []

for chain in struc.chains:

for residue in chain.residues[1:-1]:

phi, psi = getPhiPsi(residue)

phiList.append(phi)

psiList.append(psi)

pyplot.scatter(phiList, psiList) # Scatter plot

pyplot.axis([−180,180,−180,180]) # Set bounds

pyplot.show()

Note  that  if  you  have  a  number  of  structural  models,  for  example  from  an  ensemble

calculated by NMR, then you may wish to calculate the average of torsion angles. There is

a function described in

Chapter 9

which does this, given that angles are a cyclic measure

(e.g.  −180°,  180°  and  540°  mean  the  same  thing)  and  a  simple  numerical  average  is  no

good.