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