Sunday, August 1, 2010

Principle axes of inertia as a mean to reorient a molecule

In structural biology, it is often helpful to reorient a molecule in a specific view for easy visualization and comparison. With no a priori knowledge of a feature, the (objective) principle axes of inertia can be used as a mean to automatically reset a structure into its most extended orientation (and two other orthogonal views). In 3DNA, rotate_mol has such a functionality which is explored by blocview to generate the characteristic, simplified molecular images adopted by the RCSB PDB and NDB.

Recently, the thread (initiated by Dr. Pascal Auffinger) at the 3DNA forum titled "rotate_mol ambiguous orientations" has led me to think the topic deeper. As always, the issue is best illustrated with an example. The file 355d_A6_T19.pdb contains the coordinates of the base pair A6-T19 of PDB entry 355d. Let xyz be a variable in Octave (Matlab) containing the coordinates of all atoms in file 355d_A6_T19.pdb, the following m-file provides the detailed procedure:
xyz = [
13.954 12.080 15.071
13.092 11.246 15.953
14.769 11.453 14.020
13.056 13.202 14.406
12.184 13.965 15.222
11.308 14.846 14.376
12.033 15.991 13.896
10.654 14.188 13.168
9.262 14.480 13.256
11.373 14.820 11.974
11.810 16.185 12.511
13.051 16.742 11.953
14.253 16.097 11.802
15.197 16.853 11.305
14.585 18.084 11.106
15.066 19.322 10.605
16.329 19.520 10.209
14.198 20.357 10.535
12.938 20.162 10.946
12.370 19.048 11.444
13.259 18.034 11.497
14.865 29.619 5.710
14.172 30.844 5.281
16.107 29.695 6.515
13.846 28.701 6.510
12.684 28.174 5.864
11.909 27.311 6.826
12.626 26.090 7.129
11.584 27.960 8.173
10.189 27.833 8.395
12.409 27.166 9.184
12.544 25.808 8.521
13.716 24.979 8.867
13.496 23.688 9.310
12.381 23.202 9.467
14.631 22.971 9.551
15.943 23.388 9.393
16.852 22.613 9.632
16.107 24.759 8.938
17.488 25.309 8.745
14.995 25.475 8.711 ]

cov_mtx = cov(xyz)

% cov_mtx =
% 3.4094 2.3860 -1.4222
% 2.3860 33.7089 -15.7175
% -1.4222 -15.7175 7.8786

[v, lambda] = eig(cov_mtx)

% v =
% 0.094830 0.992835 -0.072709
% 0.419520 -0.106092 -0.901525
% 0.902779 -0.054989 0.426575

% lambda =
% 0.42535 0 0
% 0 3.23319 0
% 0 0 41.33836

x = v(:, 2);
y = v(:, 3);
z = cross(x, y);

% z =
% -0.094830
% -0.419520
% -0.902779

num = size(xyz, 1);
xyz2 = (xyz - ones(num, 1) * mean(xyz)) * [x y z]

% xyz2 =
% 1.148136 10.032849 -0.514731
% 0.332293 11.223635 -0.879359
% 2.081609 10.090517 0.619843
% .................................
% 2.277813 -4.170312 -0.501254
% 3.601180 -4.848890 -0.688714
% 1.110301 -4.831784 -0.491249
Among other subtitles, note that:
  1. lambda 3 has the largest value (41.33836), i.e., with the biggest variance, and its corresponding eigenvector is set as the vertical y-axis.
  2. lambda 2 shows intermediate variation (3.23319), and its corresponding eigenvector is set as the horizontal x-axis.
  3. The z-axis is derived as the cross-product of x- and y-axes to ensure a right-handed frame. Note that it is the opposite of eigenvector 1.
  4. The original xyz coordinates are moved to its geometric center [xyz - ones(num, 1) * mean(xyz)], and then projected onto the preferred frame to get the transformed coordinates (xyz2) with the molecule in the most extended view (vertical then horizontal).

1 comment:

  1. I used this code at work today, very helpful!! Thank you, your blog is so interesting!

    ReplyDelete

You are welcome to make a comment. Just remember to be specific and follow common-sense etiquette.