Drawing Ramachandran (phi/psi) plots

So you have a PDB file and you want to use Python to draw a Ramachandran diagram? This page follows on from Calculating the angles. If you look at the output from other tools we ideally want to colour the image background to show the expected regions - this makes "odd" residues really stand out.

Its very simple to write the angles to a file, and use another tool to draw a scatter plot for you - maybe with your favourite spreadsheet program. Or, if you want the background filled in, you might generate an input file for Peter Robinson's Java Program For Drawing Ramachandran Plots.

Using suitable python libraries, its probably very easy to draw a (ϕ,ψ) scatter plot - but ideally we want to highlight the "allowed" and "favoured" regions in the background of the figure. This requires a definition of those regions in a format that python can understand... and a way of using this for a contour plot.

Loading the Expected Occupancy Grid

I tried to repeat the analysis by Lovell et al (2003) (see references) of their Top 500 high resolution PDB files to build the required density profile. In the end, I was able to get in touch with the group and they sent me the relevant data (files top500-angles/pct/rama/rama500-*.data, see downloads).

The Lovell et al. (2003) rama500-*.data files store the density profiles or occupancies as a grid (or matrix) with the angles split into buckets from [-180,-178] to [178, 180] which we can represent as a grid of z values at the points ϕ and ψ ranging from -179, -177, ..., to 179.

For example, their file rama500-general.data looks like this: 

# Table name/description: "Top500 General case (not Gly, Pro, or pre-Pro) B<30"
# Number of dimensions: 2
# For each dimension, 1 to 2: lower_bound  upper_bound  number_of_bins  wrapping
#   x1: -180.0 180.0 180 true
#   x2: -180.0 180.0 180 true
# List of table coordinates and values. (Value is last number on each line.)
-179.0 -179.0 0.00782923406455425
-179.0 -177.0 0.00641357067237856
-179.0 -175.0 0.005231799492823202
-179.0 -173.0 0.004234680060073368
...
179.0 179.0 0.006881355097619223

I wrote a simple python function to load these files into Numeric arrays of floats (see downloads for the full script).

Drawing the graph

I was able to draw this in R using my own R function which uses their standard scatter routines combined with the filled.contour function. Now, its possible to use RPy (R from python) to create this function and call it from within python.

The important "trick" is this step, where after importing rpy and loading the required MASS library we get R to evaluate a command that creates an R function called ramachandran.plot - note that in python you can use triple quotes to insert a string with line-breaks: 

from rpy import *
r.library("MASS")

print "Creating R function",
r("""
ramachandran.plot <- function(x.scatter, y.scatter,
    # ...
    # code removed for brevity
    # ...
    if (frame.plot)
        box()
    if (missing(plot.title))
        title(...)
    else plot.title
    invisible()
}
""")
print "Done"

Once we have done this, we can call the new R function just like any other R command (remembering of course, that full stops in function/variable/argument names in R are translated by RPy into names with underscores in python): 

r.ramachandran_plot(x_scatter=...)

And this is the end result (see downloads for the full script):

[Ramachandran Diagram for PDB code 1HMP]

Downloads

You might want to look at the following files:

This code was tested on both Windows XP (with Python 2.3.3, R 2.3.1 and RPy 0.99.2+), and on Linux (with Python 2.4, R 2.1.1 and RPy 0.99.2).