Source code

001/*
002 *                    BioJava development code
003 *
004 * This code may be freely distributed and modified under the
005 * terms of the GNU Lesser General Public Licence.  This should
006 * be distributed with the code.  If you do not have a copy,
007 * see:
008 *
009 *      http://www.gnu.org/copyleft/lesser.html
010 *
011 * Copyright for this code is held jointly by the individual
012 * authors.  These should be listed in @author doc comments.
013 *
014 * For more information on the BioJava project and its aims,
015 * or to join the biojava-l mailing list, visit the home page
016 * at:
017 *
018 *      http://www.biojava.org/
019 *
020 */
021
022package org.biojava.nbio.structure.symmetry.geometry;
023
024import org.biojava.nbio.structure.jama.EigenvalueDecomposition;
025import org.biojava.nbio.structure.jama.Matrix;
026
027import javax.vecmath.*;
028
029/**
030 *
031 * @author Peter
032 */
033public final class SuperPosition {
034
035        public static Matrix4d superpose(Point3d[] x, Point3d[] y) {
036                //superpose x onto y
037                Point3d[] ref = clonePoint3dArray(y);
038
039                Point3d ytrans = centroid(ref);
040                ytrans.negate();
041                translate(ytrans, ref);
042
043                center(x);
044
045                // calculate quaternion from relative orientation
046                Quat4d q = quaternionOrientation(x, ref);
047
048                Matrix4d rotTrans = new Matrix4d();
049                rotTrans.set(q);
050
051                // set translational component of transformation matrix
052                ytrans.negate();
053                rotTrans.setTranslation(new Vector3d(ytrans));
054
055                // tranform coordinates
056                transform(rotTrans, x);
057
058                // TODO should include translation into transformation matrix
059
060                return rotTrans;
061        }
062
063        public static Matrix4d superposeWithTranslation(Point3d[] x, Point3d[] y) {
064                //superpose x onto y
065
066                // translate to origin
067                Point3d[] xref = clonePoint3dArray(x);
068                Point3d xtrans = centroid(xref);
069                xtrans.negate();
070                translate(xtrans, xref);
071
072                Point3d[] yref = clonePoint3dArray(y);
073                Point3d ytrans = centroid(yref);
074                ytrans.negate();
075                translate(ytrans, yref);
076
077                // calculate rotational component (rotation around origin)
078                Quat4d q = quaternionOrientation(xref, yref);
079                Matrix4d rotTrans = new Matrix4d();
080                rotTrans.set(q);
081
082                // combine with x -> origin translation
083                Matrix4d trans = new Matrix4d();
084                trans.setIdentity();
085                trans.setTranslation(new Vector3d(xtrans));
086                rotTrans.mul(rotTrans, trans);
087
088                // combine with origin -> y translation
089                ytrans.negate();
090                Matrix4d transInverse = new Matrix4d();
091                transInverse.setIdentity();
092                transInverse.setTranslation(new Vector3d(ytrans));
093                rotTrans.mul(transInverse, rotTrans);
094
095                // transform x coordinates onto y coordinate frame
096                transform(rotTrans, x);
097
098                return rotTrans;
099        }
100
101        public static Matrix4d superposeAtOrigin(Point3d[] x, Point3d[] y) {
102                Quat4d q = quaternionOrientation(x, y);
103
104                Matrix4d rotTrans = new Matrix4d();
105                rotTrans.set(q);
106
107                return rotTrans;
108        }
109
110        public static Matrix4d superposeAtOrigin(Point3d[] x, Point3d[] y, AxisAngle4d axisAngle) {
111                Quat4d q = quaternionOrientation(x, y);
112                Matrix4d rotTrans = new Matrix4d();
113                rotTrans.setIdentity();
114                rotTrans.set(q);
115                axisAngle.set(q);
116                Vector3d axis = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
117                if (axis.lengthSquared() < 1.0E-6) {
118//              System.err.println("Error: SuperPosition.superposeAtOrigin: axis vector undefined!");
119                        axisAngle.x = 0;
120                        axisAngle.y = 0;
121                        axisAngle.z = 1;
122                        axisAngle.angle = 0;
123                } else {
124                        axis.normalize();
125                        axisAngle.x = axis.x;
126                        axisAngle.y = axis.y;
127                        axisAngle.z = axis.z;
128                }
129                transform(rotTrans, x);
130//              System.out.println("SuperPosition.superposeAtOrigin: rotTrans");
131//              System.out.println(rotTrans);
132//        Matrix4d temp = new Matrix4d();
133//        temp.setIdentity();
134//        temp.set(axisAngle);
135//        System.out.println("SuperPosition.superposeAtOrigin: from axisAngle");
136//              System.out.println(temp);
137                return rotTrans;
138        }
139
140
141
142
143        public static double rmsd(Point3d[] x, Point3d[] y) {
144                double sum = 0.0;
145                for (int i = 0; i < x.length; i++) {
146                        sum += x[i].distanceSquared(y[i]);
147                }
148                return Math.sqrt(sum/x.length);
149        }
150
151        public static double rmsdMin(Point3d[] x, Point3d[] y) {
152                double sum = 0.0;
153                for (int i = 0; i < x.length; i++) {
154                        double minDist = Double.MAX_VALUE;
155                        for (int j = 0; j < y.length; j++) {
156                           minDist = Math.min(minDist, x[i].distanceSquared(y[j]));
157                        }
158                        sum += minDist;
159                }
160                return Math.sqrt(sum/x.length);
161        }
162
163        /**
164         * Returns the TM-Score for two superimposed sets of coordinates
165         * Yang Zhang and Jeffrey Skolnick, PROTEINS: Structure, Function, and Bioinformatics 57:702–710 (2004)
166         * @param x coordinate set 1
167         * @param y coordinate set 2
168         * @param lengthNative total length of native sequence
169         * @return
170         */
171        public static double TMScore(Point3d[] x, Point3d[] y, int lengthNative) {
172                double d0 = 1.24 * Math.cbrt(x.length - 15.0) - 1.8;
173                double d0Sq = d0*d0;
174
175                double sum = 0;
176                for(int i = 0; i < x.length; i++) {
177                        sum += 1.0/(1.0 + x[i].distanceSquared(y[i])/d0Sq);
178                }
179
180                return sum/lengthNative;
181        }
182
183        public static double GTSlikeScore(Point3d[] x, Point3d[] y) {
184                int contacts = 0;
185
186                for (Point3d px: x) {
187                        double minDist = Double.MAX_VALUE;
188
189                        for (Point3d py: y) {
190                           minDist = Math.min(minDist, px.distanceSquared(py));
191                        }
192
193                        if (minDist > 64) continue;
194                        contacts++;
195
196                        if (minDist > 16) continue;
197                        contacts++;
198
199                        if (minDist > 4) continue;
200                        contacts++;
201
202                        if (minDist > 1) continue;
203                        contacts++;
204                }
205
206                return contacts*25.0/x.length;
207        }
208
209        public static int contacts(Point3d[] x, Point3d[] y, double maxDistance) {
210                int contacts = 0;
211                for (int i = 0; i < x.length; i++) {
212                        double minDist = Double.MAX_VALUE;
213                        for (int j = 0; j < y.length; j++) {
214                           minDist = Math.min(minDist, x[i].distanceSquared(y[j]));
215                        }
216                        if (minDist < maxDistance*maxDistance) {
217                                contacts++;
218                        }
219                }
220                return contacts;
221        }
222
223        public static void transform(Matrix4d rotTrans, Point3d[] x) {
224                for (Point3d p: x) {
225                        rotTrans.transform(p);
226                }
227        }
228
229        public static void translate(Point3d trans, Point3d[] x) {
230                for (Point3d p: x) {
231                        p.add(trans);
232                }
233        }
234
235        public static void center(Point3d[] x) {
236                Point3d center = centroid(x);
237                center.negate();
238                translate(center, x);
239        }
240
241        public static Point3d centroid(Point3d[] x) {
242                Point3d center = new Point3d();
243                for (Point3d p: x) {
244                        center.add(p);
245                }
246                center.scale(1.0/x.length);
247                return center;
248        }
249
250        private static Quat4d quaternionOrientation(Point3d[] a, Point3d[] b)  {
251                Matrix m = calcFormMatrix(a, b);
252                EigenvalueDecomposition eig = m.eig();
253                double[][] v = eig.getV().getArray();
254                Quat4d q = new Quat4d(v[1][3], v[2][3], v[3][3], v[0][3]);
255                q.normalize();
256                q.conjugate();
257//       System.out.println("SuperPosition.quaternionOrientation: " + q);
258                return q;
259        }
260
261        private static Matrix calcFormMatrix(Point3d[] a, Point3d[] b) {
262                double xx=0.0, xy=0.0, xz=0.0, yx=0.0, yy=0.0, yz=0.0, zx=0.0, zy=0.0, zz=0.0;
263
264                for (int i = 0; i < a.length; i++) {
265                        xx += a[i].x * b[i].x;
266                        xy += a[i].x * b[i].y;
267                        xz += a[i].x * b[i].z;
268                        yx += a[i].y * b[i].x;
269                        yy += a[i].y * b[i].y;
270                        yz += a[i].y * b[i].z;
271                        zx += a[i].z * b[i].x;
272                        zy += a[i].z * b[i].y;
273                        zz += a[i].z * b[i].z;
274                }
275
276                double[][] f = new double[4][4];
277                f[0][0] = xx + yy + zz;
278                f[0][1] = zy - yz;
279                f[1][0] = f[0][1];
280                f[1][1] = xx - yy - zz;
281                f[0][2] = xz - zx;
282                f[2][0] = f[0][2];
283                f[1][2] = xy + yx;
284                f[2][1] = f[1][2];
285                f[2][2] = yy - zz - xx;
286                f[0][3] = yx - xy;
287                f[3][0] = f[0][3];
288                f[1][3] = zx + xz;
289                f[3][1] = f[1][3];
290                f[2][3] = yz + zy;
291                f[3][2] = f[2][3];
292                f[3][3] = zz - xx - yy;
293
294                return new Matrix(f);
295        }
296
297        public static Point3d[] clonePoint3dArray(Point3d[] x) {
298                Point3d[] clone = new Point3d[x.length];
299                for (int i = 0; i < x.length; i++) {
300                   clone[i] = new Point3d(x[i]);
301                }
302                return clone;
303        }
304}