How to calculate torsion angle?

Given the x-, y-, and z-coordinates of four points in 3-dimensional (3D) space, how to calculate torsion angle? Overall, this is a well-solved problem in structural biology, and one can find detailed description of the algorithm in text books and on-line documents. The algorithm for calculating torsion angle is implementated in virtually every software package in structural chemistry or biology.

Basic as it is, however, in my experience, it is very important to have a detailed appreciation of how the method works in order to really get into the 3D world. Here is a worked example using Octave/Matlab of my simplified, geometry-based implementation of how to calculate torsion angle, including how to determine its sign. No theory or (complicated) mathematical formula, just a step-by-step illustration of how I solve this problem.

1. Coordinates of four points in variable abcd:
   

   abcd = [ 21.350  31.325  22.681
   22.409  31.286  21.483
   22.840  29.751  21.498
   23.543  29.175  22.594 ];

2. Two auxiliary functions: norm_vec() to normalize a vector; get_orth_norm_vec() to get the orthogonal component (normalized) of a vector with reference to another vector, which should have been normalized.
   

   function ovec = norm_vec(vec)
     ovec = vec / norm(vec);
   endfunction
   
   function ovec = get_orth_norm_vec(vec, vref)
     temp = vec - vref * dot(vec, vref);
     ovec = norm_vec(temp);
   endfunction

3. Get three vectors: b_c is the normalized vector b→c; b_a_orth is the orthogonal component (normalized) of vector b→a with reference to b→c; c_d_orth is similarly defined, as the orthogonal component (normalized) of vector c→d with reference to b→c.
   

   b_c = norm_vec(abcd(3, :) - abcd(2, :))
     % [0.2703158  -0.9627257   0.0094077]
   b_a_orth = get_orth_norm_vec(abcd(1, :) - abcd(2, :), b_c)
     % [-0.62126  -0.16696   0.76561]
   c_d_orth = get_orth_norm_vec(abcd(4, :) - abcd(3, :), b_c)
     % [0.41330   0.12486   0.90199]

4. Now the torsion angle is the angle between the two vectors, b_a_orth and c_d_orth, and can be easily calculated by their dot product. The sign of the torsion angle is determined by the relative orientation of the cross product of the same two vectors with reference to the middle vector b→c. Here they are in opposite direction, thus the torsion angle is negative.
   

   angle_deg = acos(dot(b_a_orth, c_d_orth)) * 180 / pi
     % 65.609
   sign = dot(cross(b_a_orth, c_d_orth), b_c)
     % -0.91075
   if (sign < 0)
     ang_deg = -angle_deg  % -65.609
   endif
   

A related concept is the so-called dihedral angle, or more generally the angle between two planes. As long as the normal vectors to the two corresponding planes are defined, the angle between them is easy to work out.

Moreover, the method to calculate twist angles of helical nucleic acid structures in SCHNAaP and 3DNA is essentially the same.

Fit a least squares plane to a set of points

In the early days when I first got into the field of DNA structures, I was interested in knowing details on how the "best mean (base-pair) plane" is actually defined. I could not find a solution from the articles I read at that time. For clarification, I sent emails to some authors whose papers mentioned the concept, still the results were less than satisfactory.

The mystery was finally solved after I read the following two articles from early days (traced initially from a comment in the source code of Dickerson's NewHelix, if I recall it correctly):

• D. M. Blow (1960). ``To Fit a Plane to a Set of Points by Least Squares.'' Acta Cryst., 13, 168.

• V. Schomaker et al. (1959). ``To Fit a Plane or a Line to a Set of Points by Least Squares.'' Acta Cryst., 12, 600-604.

As is clear from the titles, they focused just on least-squares plane/line fitting algorithms. More specifically, the worked example(s) help make the underlying concept clear.

As a concrete example, let's work out the least-square plane of an adenine using Octave/Matlab. The implementation here is much simpler than those in the literature cited above. Note: this example is mentioned in the 3DNA home page in the Technical Details section, so the order of atoms is kept the same as there.

The x-, y-, and z-coordinates are as follows:

xyz = [
16.461  17.015  14.676  %  91  N9
15.775  18.188  14.459  % 100  C4
14.489  18.449  14.756  %  99  N3
14.171  19.699  14.406  %  98  C2
14.933  20.644  13.839  %  97  N1
16.223  20.352  13.555  %  95  C6
16.984  21.297  12.994  %  96  N6
16.683  19.056  13.875  %  94  C5
17.918  18.439  13.718  %  93  N7
17.734  17.239  14.207  %  92  C8
];

Cmtx = cov(xyz)
[V, LAMBDA] = eig(Cmtx)

The covariance matrix (Cmtx), its corresponding eigenvectors (V) and eigenvalues (LAMBDA) are as follows:

Cmtx =
1.66800  -0.50154  -0.32530
-0.50154   2.06703  -0.58401
-0.32530  -0.58401   0.30605

V =
0.27367   0.83264  -0.48147
0.32243   0.39220   0.86152
0.90617  -0.39102  -0.16113

LAMBDA =
0.00001   0.00000   0.00000
0.00000   1.58452   0.00000
0.00000   0.00000   2.45656

The smallest eigenvalue of Cmtx, 0.00001 is close to zeros, indicating that the base is almost perfectly planar. The corresponding unit eigenvector, i.e., the least squares base normal is: [0.27367 0.32243 0.90617]. Of course, the geometric average of the atoms defines a point the ls-plane pass through.

Referring back to the concept of "best mean (base-pair) plane", it could simply be defined using all base atoms of the pair, or alternatively, as the mean of the two base normals. Obviously, they would lead to (slight) numerical differences. Astute viewers would also notice some other subtleties.

Eddy's primer on the BLOSUM62 alignment score matrix

As a Linux/Unix fun, I like its philosophy, as summarized by Doug McIlroy, very much: "Write programs that do one thing and do it well. Write programs to work together." In science, I enjoy more reading an article that focuses on one point and address it clearly and thoroughly. It is only after a complete understanding of the components can it be possible to combine them in unique, purpose-specific ways. Those days, though, such type of simple-and-clean articles is no longer that common as it was in the early days, say 1960s or 70s. This post is the first of a series on such articles I've found useful, or on computational tricks that I have learned over the years.

I came across the primer titled "Where did the BLOSUM62 alignment score matrix come from?" by Sean Eddy [Nat Biotechnol. 2004 Aug;22(8):1035-6] early this month through the BioConductor mailing list where it was recommended by Dr. Philipp Pagel as a "well written article". I read the title and the short abstract from PubMed, and then download the PDF version of the whole article.

As a primer, it is only two pages long (short), and reading twice won't take that much time. It explains the meaning of the BLOSUM62 amino acid score matrix and where is comes from clearly. The number 62, for example, stands for a threshold of 62% identity. Other percentages, such as 80% (more Conservative), 45% (more divergent) are also possible and may be more suitable for specific applications. As noted by Eddy, "Empirically, the BLOSUM matrices have performed very well. BLOSUM62 has become a de facto standard for many protein alignment programs."

With the above background, it is much easier to understand the fundamental difference of the default scoring systems between FASTA/WU-BLASTN vs NCBI BLASTN for alignment of DNA sequences: the former is optimal for alignments at the ‘twilight zone’ (65% identity), while the later (NCBI BLASTN) is optimal for homologous DNA at a much higher 95% identity level. Knowing of such subtlety is very important in avoid making false conclusions.

More significantly (to me), there is a supplemental material -- a well-documented, self-contained "ANSI C program for calculating the implicit target frequencies p_ab of a score matrix": it clarifies every details for those who want to get to the bottom of the topic.