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.core;
023
024import org.biojava.nbio.structure.geometry.SuperPositions;
025
026import javax.vecmath.Matrix4d;
027import javax.vecmath.Point3d;
028
029import java.util.List;
030
031/**
032 *
033 * @author Peter
034 */
035public class QuatSuperpositionScorer {
036
037        /**
038         * Returns minimum, mean, and maximum RMSD and TM-Score for two superimposed sets of subunits
039         *
040         * TM score: Yang Zhang and Jeffrey Skolnick, PROTEINS: Structure, Function, and Bioinformatics 57:702–710 (2004)
041         * @param subunits subunits to be scored
042         * @param transformation transformation matrix
043         * @param permutations permutation that determines which subunits are superposed
044         * @return
045         */
046        public static QuatSymmetryScores calcScores(QuatSymmetrySubunits subunits, Matrix4d transformation, List<Integer> permutation) {
047                QuatSymmetryScores scores = new QuatSymmetryScores();
048
049                double minTm = Double.MAX_VALUE;
050                double maxTm = Double.MIN_VALUE;
051                double minRmsd = Double.MAX_VALUE;
052                double maxRmsd = Double.MIN_VALUE;
053
054                double totalSumTm = 0;
055                double totalSumDsq = 0;
056                double totalLength = 0;
057
058                Point3d t = new Point3d();
059                List<Point3d[]> traces = subunits.getTraces();
060
061                // loop over the Calpha atoms of all subunits
062                for (int i = 0; i < traces.size(); i++) {
063                        // in helical systems not all permutations involve all subunit. -1 indicates subunits that should not be permuted.
064                         if (permutation.get(i) == -1) {
065                                continue;
066                         }
067                        // get original subunit
068                        Point3d[] orig = traces.get(i);
069                        totalLength += orig.length;
070
071                        // get permuted subunit
072                        Point3d[] perm = traces.get(permutation.get(i));
073
074                        // calculate TM specific parameters
075                        int tmLen = Math.max(orig.length, 17);  // don't let d0 get negative with short sequences
076                        double d0 = 1.24 * Math.cbrt(tmLen - 15.0) - 1.8;
077                        double d0Sq = d0 * d0;
078
079                        double sumTm = 0;
080                        double sumDsq = 0;
081                        for (int j = 0; j < orig.length; j++) {
082                                // transform coordinates of the permuted subunit
083                                t.set(perm[j]);
084                                transformation.transform(t);
085
086                                double dSq = orig[j].distanceSquared(t);
087                                sumTm += 1.0/(1.0 + dSq/d0Sq);
088                                sumDsq += dSq;
089                        }
090
091                        // scores for individual subunits
092                        double sTm = sumTm/tmLen;
093                        minTm = Math.min(minTm, sTm);
094                        maxTm = Math.max(maxTm, sTm);
095
096                        double sRmsd = Math.sqrt(sumDsq/orig.length);
097                        minRmsd = Math.min(minRmsd,  sRmsd);
098                        maxRmsd = Math.max(maxRmsd,  sRmsd);
099
100                        totalSumTm += sumTm;
101                        totalSumDsq += sumDsq;
102                }
103
104                // save scores for individual subunits
105                scores.setMinRmsd(minRmsd);
106                scores.setMaxRmsd(maxRmsd);
107                scores.setMinTm(minTm);
108                scores.setMaxTm(maxTm);
109
110                // save mean scores over all subunits
111                scores.setTm(totalSumTm/totalLength);
112                scores.setRmsd(Math.sqrt(totalSumDsq/totalLength));
113
114                // add intra subunit scores
115                calcIntrasubunitScores(subunits, transformation, permutation, scores);
116
117                return scores;
118        }
119
120        private static void calcIntrasubunitScores(QuatSymmetrySubunits subunits, Matrix4d transformation, List<Integer> permutation, QuatSymmetryScores scores) {
121                double totalSumTm = 0;
122                double totalSumDsq = 0;
123                double totalLength = 0;
124
125                List<Point3d[]> traces = subunits.getTraces();
126
127                // loop over the Calpha atoms of all subunits
128                for (int i = 0; i < traces.size(); i++) {
129                        // in helical systems not all permutations involve all subunit. -1 indicates subunits that should not be permuted.
130                         if (permutation.get(i) == -1) {
131                                continue;
132                         }
133                        // get original subunit
134                        Point3d[] orig = traces.get(i);
135                        totalLength += orig.length;
136
137                        // get permuted subunit
138                        Point3d[] perm = traces.get(permutation.get(i));
139
140                        // calculate TM specific parameters
141                        int tmLen = Math.max(orig.length, 17);  // don't let d0 get negative with short sequences
142                        double d0 = 1.24 * Math.cbrt(tmLen - 15.0) - 1.8;
143                        double d0Sq = d0 * d0;
144
145                        double sumTm = 0;
146                        double sumDsq = 0;
147                        Point3d[] trans = new Point3d[orig.length];
148                        for (int j = 0; j < orig.length; j++) {
149                                trans[j] = new Point3d(perm[j]);
150                        }
151
152                        // superpose individual subunits
153                        SuperPositions.superposeAndTransform(orig, trans);
154                        for (int j = 0; j < orig.length; j++) {
155                                double dSq = orig[j].distanceSquared(trans[j]);
156                           sumTm += 1.0/(1.0 + dSq/d0Sq);
157                           sumDsq += dSq;
158                        }
159
160                        totalSumTm += sumTm;
161                        totalSumDsq += sumDsq;
162                }
163                // save mean scores over all subunits
164                scores.setRmsdIntra(Math.sqrt(totalSumDsq/totalLength));
165                scores.setTmIntra(totalSumTm/totalLength);
166        }
167
168}