Kabsch algorithm
The Kabsch algorithm, also known as the Kabsch-Umeyama algorithm,[1] named after Wolfgang Kabsch and Shinji Umeyama, is a method for calculating the optimal rotation matrix that minimizes the RMSD (root mean squared deviation) between two paired sets of points. It is useful for point-set registration in computer graphics, and in cheminformatics and bioinformatics to compare molecular and protein structures (in particular, see root-mean-square deviation (bioinformatics)). The algorithm only computes the rotation matrix, but it also requires the computation of a translation vector. When both the translation and rotation are actually performed, the algorithm is sometimes called partial Procrustes superimposition (see also orthogonal Procrustes problem).
Description
Let P and Q be two sets, each containing N points in [math]\displaystyle{ \mathbb{R}^n }[/math]. For simplicity, we will consider the three-dimensional case ([math]\displaystyle{ n = 3 }[/math]). The sets P and Q can each be represented by N × 3 matrices with the first row containing the coordinates of the first point, the second row containing the coordinates of the second point, and so on, as shown in this matrix:
[math]\displaystyle{ \begin{pmatrix} x_1 & y_1 & z_1 \\ x_2 & y_2 & z_2 \\ \vdots & \vdots & \vdots \\ x_N & y_N & z_N \end{pmatrix} }[/math]
The algorithm works in three steps: a translation, the computation of a covariance matrix, and the computation of the optimal rotation matrix.
Translation
Both sets of coordinates must be translated first, so that their centroid coincides with the origin of the coordinate system. This is done by subtracting from the point coordinates of the respective centroid.
Computation of the covariance matrix
The second step consists of calculating a matrix H. In matrix notation,
- [math]\displaystyle{ H = P^\mathsf{T}Q \, }[/math]
or, using summation notation,
- [math]\displaystyle{ H_{ij} = \sum_{k = 1}^N P_{ki} Q_{kj}, }[/math]
which is a cross-covariance matrix when P and Q are seen as data matrices.
Computation of the optimal rotation matrix
It is possible to calculate the optimal rotation R based on the matrix formula
- [math]\displaystyle{ R = \left(H^\mathsf{T} H\right)^\frac12 H^{-1} }[/math]
but implementing a numerical solution to this formula becomes complicated when all special cases are accounted for (for example, the case of H not having an inverse).
If singular value decomposition (SVD) routines are available, the optimal rotation, R, can be calculated using the following simple algorithm.
First, calculate the SVD of the covariance matrix H.
- [math]\displaystyle{ H = U \Sigma V^\mathsf{T} }[/math]
Next, decide whether we need to correct our rotation matrix to ensure a right-handed coordinate system
- [math]\displaystyle{ d = \mathrm{sign}\left(\det\left(V U^\mathsf{T}\right)\right) }[/math]
Finally, calculate our optimal rotation matrix, R, as
- [math]\displaystyle{ R = V \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & d \end{pmatrix} U^\mathsf{T} }[/math]
The optimal rotation matrix can also be expressed in terms of quaternions.[2][3][4][5] This alternative description has been used in the development of a rigorous method for removing rigid-body motions from molecular dynamics trajectories of flexible molecules.[6] In 2002 a generalization for the application to probability distributions (continuous or not) was also proposed.[7]
Generalizations
The algorithm was described for points in a three-dimensional space. The generalization to D dimensions is immediate.
External links
This SVD algorithm is described in more detail at https://web.archive.org/web/20140225050055/http://cnx.org/content/m11608/latest/
A Matlab function is available at http://www.mathworks.com/matlabcentral/fileexchange/25746-kabsch-algorithm
A C++ implementation (and unit test) using Eigen
A Python script is available at https://github.com/charnley/rmsd. Another implementation can be found in SciPy.
A free PyMol plugin easily implementing Kabsch is [1]. (This previously linked to CEalign [2], but this uses the Combinatorial Extension (CE) algorithm.) VMD uses the Kabsch algorithm for its alignment.
The FoldX modeling toolsuite incorporates the Kabsch algorithm to measure RMSD between Wild Type and Mutated protein structures.
See also
References
- ↑ Lawrence, Jim; Bernal, Javier; Witzgall, Christoph (2019-10-09). "A Purely Algebraic Justification of the Kabsch-Umeyama Algorithm" (in en). Journal of Research of the National Institute of Standards and Technology 124: 124028. doi:10.6028/jres.124.028. ISSN 2165-7254. PMID 34877177. PMC 7340555. https://nvlpubs.nist.gov/nistpubs/jres/124/jres.124.028.pdf.
- ↑ Horn, Berthold K. P. (1987-04-01). "Closed-form solution of absolute orientation using unit quaternions" (in EN). Journal of the Optical Society of America A 4 (4): 629. doi:10.1364/josaa.4.000629. ISSN 1520-8532. Bibcode: 1987JOSAA...4..629H.
- ↑ Kneller, Gerald R. (1991-05-01). "Superposition of Molecular Structures using Quaternions". Molecular Simulation 7 (1–2): 113–119. doi:10.1080/08927029108022453. ISSN 0892-7022.
- ↑ Coutsias, E. A.; Seok, C.; Dill, K. A. (2004). "Using quaternions to calculate RMSD". J. Comput. Chem. 25 (15): 1849–1857. doi:10.1002/jcc.20110. PMID 15376254.
- ↑ Petitjean, M. (1999). "On the root mean square quantitative chirality and quantitative symmetry measures". J. Math. Phys. 40 (9): 4587–4595. doi:10.1063/1.532988. Bibcode: 1999JMP....40.4587P. https://hal.archives-ouvertes.fr/hal-02122820/file/PMP.JMP_1999.pdf.
- ↑ Chevrot, Guillaume; Calligari, Paolo; Hinsen, Konrad; Kneller, Gerald R. (2011-08-24). "Least constraint approach to the extraction of internal motions from molecular dynamics trajectories of flexible macromolecules". J. Chem. Phys. 135 (8): 084110. doi:10.1063/1.3626275. ISSN 0021-9606. PMID 21895162. Bibcode: 2011JChPh.135h4110C.
- ↑ Petitjean, M. (2002). "Chiral mixtures". J. Math. Phys. 43 (8): 4147–4157. doi:10.1063/1.1484559. Bibcode: 2002JMP....43.4147P. https://hal.archives-ouvertes.fr/hal-02122882/file/PMP.JMP_2002.pdf.
- Kabsch, Wolfgang (1976). "A solution for the best rotation to relate two sets of vectors". Acta Crystallographica A32 (5): 922. doi:10.1107/S0567739476001873. Bibcode: 1976AcCrA..32..922K.
- With a correction in Kabsch, Wolfgang (1978). "A discussion of the solution for the best rotation to relate two sets of vectors". Acta Crystallographica A34 (5): 827–828. doi:10.1107/S0567739478001680. Bibcode: 1978AcCrA..34..827K.
- Lin, Ying-Hung; Chang, Hsun-Chang; Lin, Yaw-Ling (December 15–17, 2004). "A Study on Tools and Algorithms for 3-D Protein Structures Alignment and Comparison". International Computer Symposium. Taipei, Taiwan.
- Umeyama, Shinji (1991). "Least-Squares Estimation of Transformation Parameters Between Two Point Patterns". IEEE Trans. Pattern Anal. Mach. Intell. 13 (4): 376–380. doi:10.1109/34.88573.
 |  |
