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.
Do'stlaringiz bilan baham: |