Hi, Pavel, Yes, my suggestion would be specific for double-stranded nucleic acid. A separate but similar one would be required for protein (I've written and used that one). For protein the bf straight line would not be so good because its a single helix, so at the end atoms are on only one side. Say you have 1.5 turns of helix. the best fit line would be skewed to go through the ends of the backbone (as Tim Gruene noted). For a double helix it would be less of a problem since you have one chain on each side of the axis to balance. But in the extreme, say when the length is less than the thickness (i.e. a disk), the best-fit line could be closer to the diameter than the axis. I wouldn't pretend to advise you on programming linear algerba, but for my own edification I will study your matrix and try to understand, with help from the page Tim linked. At first glance it looks like S.V. decomposition on n vectors in 3 space, or 3 vectors in N space, where the singular vectors, values are taken as the eigenvectors, values of the matrix U'V or UV'. Or the normal equations for solving a least-squares problem with a generalized inverse involving [M'M]^-1 but then you would go on to invert the matrix rather than finding its eigenvectors. Pavel Afonine wrote:
Hi Ed,
interesting idea! Although I was thinking to have a tool that is a little more general and a little less context dependent. Say you have two clouds of points that are (thinking in terms of macromolecules) two alpha helices (for instance), and you want to know the angle between the axes of the two helices. How would I approach this?..
First, for each helix I would compute a symmetric 3x3 matrix like this:
sum(xn-xc)**2 sum(xn-xc)*(yn-xc) sum(xn-xc)*(zn-zc) sum(xn-xc)*(yn-xc) sum(yn-yc)**2 sum(yn-yc)*(yz-zc) sum(xn-xc)*(zn-zc) sum(yn-yc)*(yz-zc) sum(zn-zc)**2
where (xn,yn,zn) is the coordinate of nth atom, the sum is taken over all atoms, and (xc,yc,zc) is the coordinate of the center of mass.
Second, for each of the two matrices I would find its eigen-values and eigen-vectors, and select eigen-vectors corresponding to largest eigenvalues.
Finally, the desired angle is the angle between the two eigen-vectors found above, which is computed trivially. I think this a little simpler than finding the best fit for a 3D line.
What you think?
Pavel
On 1/20/14, 2:14 PM, Edward A. Berry wrote:
Pavel Afonine wrote: . .
The underlying procedure would do the following: - extract two sets of coordinates of atoms corresponding to two provided atom selections; - draw two optimal lines (LS fit) passing through the above sets of coordinates; - compute and report angle between those two lines?
This could be innacurate for very short helices (admittedly not the case one usually would be looking for angles), or determining the axis of a short portion of a curved helix. A more accurate way to determine the axis- have a long canonical duplex constructed with its axis along Z (0,0,1). Superimpose as many residues of that as required on the duplex being tested, using only backbone atoms or even only phosphates. Operate on (0,0,1) with the resulting operator (i.e. take the third column of the rotation matrix) and use that as a vector parallel to the axis of the duplex being tested. _______________________________________________ phenixbb mailing list [email protected] http://phenix-online.org/mailman/listinfo/phenixbb
_______________________________________________ phenixbb mailing list [email protected] http://phenix-online.org/mailman/listinfo/phenixbb