protein_analysis.py

[Konrad Hinsen]

# This example demonstrates some of the analysis features of MMTK.
# Two protein structures from the PDB are compared in some detail.

# Disclaimer: I know nothing about this protein, and the analysis
# I present here may be of no biological interest.

from mmtk import *

# First we read the two PDB files.

sequences_1 = PDBFile('4q21.pdb').readSequenceWithConfiguration()
sequences_2 = PDBFile('6q21.pdb').readSequenceWithConfiguration()

# The first file contains a monomer, the second one a tetramer in
# which each chain is almost identical to the monomer from the first
# file. We have to cut off the last (incomplete) residue from the
# monomer and the last three residues of each chain of the tetramers
# to get matching sequences. We'll just deal with one of the chains of
# the tetramer here.

monomer = sequences_1[0][:-1]
tetramer_1 = sequences_2[0][:-3]

# Next we build a PeptideChain for the monomer. We have to specify
# that we have no C terminus (missing) and no hydrogens.

chain = PeptideChain(monomer, hydrogens=0, c_terminus=0)

# How big is this beast? For a quick guess, print the size of the smallest
# rectangular box containing it:

p1, p2 = chain.boundingBox()
print "Size of a bounding box: %4.2f x %4.2f x %4.2f nm" % tuple(p2-p1)

# Formalities: we define a universe and put the chain inside.
# Then we make a copy of the configuration of the universe for later use.
# This is necessary because configurations are defined only for
# universes (for various good reasons...).

universe = InfiniteUniverse()
universe.addObject(chain)
monomer_configuration = copy(universe.configuration())

# Next we change the conformation of the protein to that of the first
# tetramer chain.

chain.setConfiguration(tetramer_1)

# Now we want to get rid of the difference in position and orientation
# of the two configurations, so we calculate the transformation that
# minimizes the RMS difference.

transformation, rms = chain.findTransformation(monomer_configuration)

# That's the transformation *from* the current *to* the monomer
# configuration. Let's analyze it:

print "Translation/rotation fit:"
print "  RMS difference:", rms

translation = transformation.translation()
rotation = transformation.rotation()
axis, angle = rotation.axisAndAngle()

print "  Rotation axis: ", axis
print "  Rotation angle: ", angle
print "  Translation: ", translation.displacement()

# Note that order is important: the rotation (around the origin) is
# applied first, and then the translation. If you want to do it the
# other way round, the translation displacement must be different.

# Note also the units: distances in nanometers, angles in radians!
# If you insist on Angstrom and degrees, you'll have to type a bit more:

print "  Rotation angle in degrees: ", angle/Units.deg
print "  Translation in Angstrom: ", translation.displacement()/Units.Ang

# Just for fun we'll do it again, but minimizing the RMS for the
# C_alpha atoms only:

c_alphas = chain.residues().map(lambda residue: residue.peptide.C_alpha)
calpha_transformation, rms = c_alphas.findTransformation(monomer_configuration)

print "C_alpha fit:"
print "  RMS difference:", rms

translation = calpha_transformation.translation()
rotation = calpha_transformation.rotation()
axis, angle = rotation.axisAndAngle()

print "  Rotation axis: ", axis
print "  Rotation angle: ", angle
print "  Translation: ", translation.displacement()

# As expected there is little difference. Now we apply the transformation
# to the current configuration to align it with the reference configuration:

chain.applyTransformation(transformation)

# Action: We show an animation of the two configurations. This requires
# an animation-capable external viewer, e.g. VMD.

viewSequence(chain, [universe.configuration(), monomer_configuration])

# Interesting: there's a floppy loop with residues Gln61, Glu62, and Glu63
# at the tip. It seems to do a funny rotation relative to a closeby helix
# (His94-Lys101) when switching between the two configurations. Let's
# find out what this rotation is!

# First, we define the parts of the chain we are interested in. We don't
# want the sidechains to mess up the picture, so we use just the backbone:

tip = chain[60:64].backbone()
helix = chain[93:101].backbone()

# Then we align the two configurations to make the helix match as well
# as possible:

transformation, rms = helix.findTransformation(monomer_configuration)
chain.applyTransformation(transformation)

# Note that we *fit* only the helix, but *apply* the transformation
# to the whole chain.

# And now another fit for the tip of the loop. This time we don't want
# to apply the transformation, but print out its characteristics:

transformation, rms = tip.findTransformation(monomer_configuration)

print "Movement of the tip of the loop:"
print "  RMS difference:", rms

translation = transformation.translation()
rotation = transformation.rotation()
axis, angle = rotation.axisAndAngle()

print "  Rotation axis: ", axis
print "  Rotation angle: ", angle
print "  Translation: ", translation.displacement()