Tuesday, 18 May 2010

What does the pairwise RMSD of all possible conformations look like?

I have been doing some work on conformer generation and trying to think about how to measure the quality of one method over another. Since I'm a great believer in looking at the data (see also this post), I decided to calculate the pairwise RMSD of the systematically generated conformers of a molecule with 3 rotatable bonds, and then visualise these data with classical multidimensional scaling (cmdscale in R):
It turned out to be quite a nice graph. Each point here represents a different conformer. Of the three torsion angles, it seems that two were incremented in steps of 30°, while for the third Δang was 120°. Each torsion angle had a different magnitude of effect on the RMSD (or at least, this is how I interpret this graph - I haven't drilled down into the figures yet); presumably the torsion angle closer to the centre of the molecule gives rise to the three clusters of values.

The points in red are the lowest energy conformers (by MMFF94). In this graph it is difficult to see what the correspondence between these conformers is, although it is clear that the overall RMSD between some of these conformers is high.

All data was generated using a locally-modified development version of OpenBabel. Here is the R code I used:
<- read.table("tmp_dist.txt")
location <- cmdscale(m)

mycol <- read.table("tmp_en.txt")
max_scaled <- max(scaled)
min_scaled <- min(scaled)
scaled <- (-min_scaled + scaled) / (max_scaled - min_scaled)

# colmap <- rainbow(11)
lookup <- function(x) {
        # colmap[floor( (x + 0.15)*10) ]
        if(x < 0.1) {"red"}
        else {"black"}
}
colors = apply(scaled, 1, lookup)
plot(location, col = colors, pch=19)

3 comments:

Evan Bolton said...

Reading this post left me with lots of questions... :-)

How exactly did you do your RMSD? Was it simple pairwise comparison? Did you allow rigid body rotation and translation to minimize the RMSD between conformers? Did you consider automorphs (symmetry equivalent atoms)?

baoilleach said...

The RMSD was after Kabsch alignment and it included hydrogen atoms (these are not usually included in 'real life'). Automorphs were not considered - at some point I may need to sort this out.

Evan Bolton said...

Good that you minimized the RMSD first. I suspect your graphs may look a little different after considering automorphs, as many hydrogen atoms are symmetry equivalent to others. You could try heavy atoms only and see if the plot changes much. :-)