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 */
021package org.biojava.nbio.structure.symmetry.core;
022
023import org.biojava.nbio.structure.symmetry.geometry.SuperPosition;
024
025import javax.vecmath.AxisAngle4d;
026import javax.vecmath.Matrix4d;
027import javax.vecmath.Point3d;
028import javax.vecmath.Vector3d;
029
030import java.util.*;
031import java.util.Map.Entry;
032
033public class HelixSolver {
034        private Subunits subunits = null;
035        private int fold = 1;
036        private HelixLayers helixLayers = new HelixLayers();
037        private QuatSymmetryParameters parameters = null;
038        boolean modified = true;
039
040        public HelixSolver(Subunits subunits, int fold, QuatSymmetryParameters parameters) {
041                this.subunits = subunits;
042                this.fold = fold;
043                this.parameters = parameters;
044        }
045
046        public HelixLayers getSymmetryOperations() {
047                if (modified) {
048                        solve();
049                        modified = false;
050                }
051                return helixLayers;
052        }
053
054        private void solve() {
055                if (! preCheck()) {
056                        return;
057                }
058
059                HelicalRepeatUnit unit = new HelicalRepeatUnit(subunits);
060                List<Point3d> repeatUnitCenters = unit.getRepeatUnitCenters();
061                List<Point3d[]> repeatUnits = unit.getRepeatUnits();
062                Set<List<Integer>> permutations = new HashSet<List<Integer>>();
063
064                double minRise = parameters.getMinimumHelixRise() * fold;  // for n-start helix, the rise must be steeper
065                Map<Integer[], Integer> interactionMap = unit.getInteractingRepeatUnits();
066
067                int maxLayerLineLength = 0;
068
069                for (Entry<Integer[], Integer> entry : interactionMap.entrySet()) {
070                        Integer[] pair = entry.getKey();
071                        if (parameters.isVerbose()) {
072                                System.out.println("HelixSolver: pair: " + Arrays.toString(pair));
073                        }
074                        int contacts = entry.getValue();
075                        Point3d[] h1 = null;
076                        Point3d[] h2 = null;
077
078                        // trial superposition of repeat unit pairs to get a seed permutation
079                        h1 = SuperPosition.clonePoint3dArray(repeatUnits.get(pair[0]));
080                        h2 = SuperPosition.clonePoint3dArray(repeatUnits.get(pair[1]));
081                        Matrix4d transformation = SuperPosition.superposeWithTranslation(h1, h2);
082
083                        double rmsd = SuperPosition.rmsd(h1, h2);
084                        double rise = getRise(transformation, repeatUnitCenters.get(pair[0]), repeatUnitCenters.get(pair[1]));
085                        double angle = getAngle(transformation);
086
087                        if (parameters.isVerbose()) {
088                                System.out.println();
089                                System.out.println("Original rmsd: " + rmsd);
090                                System.out.println("Original rise: " + rise);
091                                System.out.println("Original angle: " + Math.toDegrees(angle));
092                        }
093
094                        if (rmsd > parameters.getRmsdThreshold()) {
095                                continue;
096                        }
097
098                        if (Math.abs(rise) < minRise) {
099                                continue;
100                        }
101
102                        // determine which subunits are permuted by the transformation
103                        List<Integer> permutation = getPermutation(transformation);
104
105                        // check permutations for validity
106
107                        // don't save redundant permutations
108                        if (permutations.contains(permutation)) {
109                                continue;
110                        }
111                        permutations.add(permutation);
112                        if (parameters.isVerbose()) {
113                                System.out.println("Permutation: " + permutation);
114                        }
115
116                        // keep track of which subunits are permuted
117                        Set<Integer> permSet = new HashSet<Integer>();
118                        int count = 0;
119                        boolean valid = true;
120                        for (int i = 0; i < permutation.size(); i++) {
121                                if (permutation.get(i) == i) {
122                                        valid = false;
123                                        break;
124                                }
125                                if (permutation.get(i) != -1) {
126                                        permSet.add(permutation.get(i));
127                                        permSet.add(i);
128                                        count++;
129                                }
130
131                        }
132
133                        // a helix a repeat unit cannot map onto itself
134                        if (! valid) {
135                                if (parameters.isVerbose()) {
136                                        System.out.println("Invalid mapping");
137                                }
138                                continue;
139                        }
140
141                        // all subunits must be involved in a permutation
142                        if (permSet.size() != subunits.getSubunitCount()) {
143                                if (parameters.isVerbose()) {
144                                        System.out.println("Not all subunits involved in permutation");
145                                }
146                                continue;
147                        }
148
149                        // if all subunit permutation values are set, then it can't be helical symmetry (must be cyclic symmetry)
150                        if (count == permutation.size()) {
151                                continue;
152                        }
153
154                        // superpose all permuted subunits
155                        List<Point3d> point1 = new ArrayList<Point3d>();
156                        List<Point3d> point2 = new ArrayList<Point3d>();
157                        List<Point3d> centers = subunits.getOriginalCenters();
158                        for (int j = 0; j < permutation.size(); j++) {
159                                if (permutation.get(j) != -1) {
160                                        point1.add(new Point3d(centers.get(j)));
161                                        point2.add(new Point3d(centers.get(permutation.get(j))));
162                                }
163                        }
164
165                        h1 = new Point3d[point1.size()];
166                        h2 = new Point3d[point2.size()];
167                        point1.toArray(h1);
168                        point2.toArray(h2);
169
170                        // calculate subunit rmsd if at least 3 subunits are available
171                        double subunitRmsd = 0;
172                        if (point1.size() > 2) {
173                                transformation = SuperPosition.superposeWithTranslation(h1, h2);
174
175                                subunitRmsd = SuperPosition.rmsd(h1, h2);
176                                rise = getRise(transformation, repeatUnitCenters.get(pair[0]), repeatUnitCenters.get(pair[1]));
177                                angle = getAngle(transformation);
178
179                                if (parameters.isVerbose()) {
180                                        System.out.println("Subunit rmsd: " + subunitRmsd);
181                                        System.out.println("Subunit rise: " + rise);
182                                        System.out.println("Subunit angle: " + Math.toDegrees(angle));
183                                }
184
185                                if (subunitRmsd > parameters.getRmsdThreshold()) {
186                                        continue;
187                                }
188
189                                if (Math.abs(rise) < minRise) {
190                                        continue;
191                                }
192
193                                if (subunitRmsd > parameters.getHelixRmsdToRiseRatio()*Math.abs(rise)) {
194                                        continue;
195                                }
196                        }
197
198                        // superpose all C alpha traces
199                        point1.clear();
200                        point2.clear();
201                        List<Point3d[]> traces = subunits.getTraces();
202                        for (int j = 0; j < permutation.size(); j++) {
203                                if (permutation.get(j) == -1) {
204                                        continue;
205                                }
206                                for (Point3d p: traces.get(j)) {
207                                        point1.add(new Point3d(p));
208                                }
209                                for (Point3d p: traces.get(permutation.get(j))) {
210                                        point2.add(new Point3d(p));
211                                }
212                        }
213
214                        h1 = new Point3d[point1.size()];
215                        h2 = new Point3d[point2.size()];
216                        point1.toArray(h1);
217                        point2.toArray(h2);
218                        Point3d[] h3 = SuperPosition.clonePoint3dArray(h1);
219                        transformation = SuperPosition.superposeWithTranslation(h1, h2);
220
221                        Point3d xtrans = SuperPosition.centroid(h3);
222                        xtrans.negate();
223
224
225                        double traceRmsd = SuperPosition.rmsd(h1, h2);
226
227                        rise = getRise(transformation, repeatUnitCenters.get(pair[0]), repeatUnitCenters.get(pair[1]));
228                        angle = getAngle(transformation);
229
230                        if (parameters.isVerbose()) {
231                                System.out.println("Trace rmsd: " + traceRmsd);
232                                System.out.println("Trace rise: " + rise);
233                                System.out.println("Trace angle: " + Math.toDegrees(angle));
234                                System.out.println("Permutation: " + permutation);
235                        }
236
237                        if (traceRmsd > parameters.getRmsdThreshold()) {
238                                continue;
239                        }
240
241                        if (Math.abs(rise) < minRise) {
242                                continue;
243                        }
244
245                        // This prevents translational repeats to be counted as helices
246                        if (angle < Math.toRadians(parameters.getMinimumHelixAngle())) {
247                                continue;
248                        }
249
250                        if (traceRmsd > parameters.getHelixRmsdToRiseRatio()*Math.abs(rise)) {
251                                continue;
252                        }
253
254                        AxisAngle4d a1 = new AxisAngle4d();
255                        a1.set(transformation);
256
257                        // save this helix rot-translation
258                        Helix helix = new Helix();
259                        helix.setTransformation(transformation);
260                        helix.setPermutation(permutation);
261                        helix.setRise(rise);
262                        // Old version of Vecmath on LINUX doesn't set element m33 to 1. Here we make sure it's 1.
263                        transformation.setElement(3, 3, 1.0);
264                        transformation.invert();
265                        QuatSymmetryScores scores = QuatSuperpositionScorer.calcScores(subunits, transformation, permutation);
266                        scores.setRmsdCenters(subunitRmsd);
267                        helix.setScores(scores);
268                        helix.setFold(fold);
269                        helix.setContacts(contacts);
270                        helix.setRepeatUnits(unit.getRepeatUnitIndices());
271                        if (parameters.isVerbose()) {
272                                System.out.println("Layerlines: " + helix.getLayerLines());
273                        }
274                        for (List<Integer> line: helix.getLayerLines()) {
275                                maxLayerLineLength = Math.max(maxLayerLineLength, line.size());
276                        }
277
278                        // TODO
279        //              checkSelfLimitingHelix(helix);
280
281                        helixLayers.addHelix(helix);
282
283                }
284                if (maxLayerLineLength < 3) {
285//                      System.out.println("maxLayerLineLength: " + maxLayerLineLength);
286                        helixLayers.clear();
287                }
288
289                return;
290        }
291
292        @SuppressWarnings("unused")
293        private void checkSelfLimitingHelix(Helix helix) {
294                HelixExtender he = new HelixExtender(subunits, helix);
295                Point3d[] extendedHelix = he.extendHelix(1);
296
297                int overlap1 = 0;
298                for (Point3d[] trace: subunits.getTraces()) {
299                        for (Point3d pt: trace) {
300                                for (Point3d pe: extendedHelix) {
301                                        if (pt.distance(pe) < 5.0) {
302                                                overlap1++;
303                                        }
304                                }
305                        }
306                }
307
308                extendedHelix = he.extendHelix(-1);
309
310                int overlap2 = 0;
311                for (Point3d[] trace: subunits.getTraces()) {
312                        for (Point3d pt: trace) {
313                                for (Point3d pe: extendedHelix) {
314                                        if (pt.distance(pe) < 3.0) {
315                                                overlap2++;
316                                        }
317                                }
318                        }
319                }
320                System.out.println("SelfLimiting helix: " + overlap1 + ", " + overlap2);
321        }
322
323        private boolean preCheck() {
324                if (subunits.getSubunitCount() < 3) {
325                        return false;
326                }
327                List<Integer> folds = this.subunits.getFolds();
328                int maxFold = folds.get(folds.size()-1);
329                return maxFold > 1;
330        }
331
332        /**
333         * Returns a permutation of subunit indices for the given helix transformation.
334         * An index of -1 is used to indicate subunits that do not superpose onto any other subunit.
335         * @param transformation
336         * @return
337         */
338        private List<Integer> getPermutation(Matrix4d transformation) {
339                double rmsdThresholdSq = Math.pow(this.parameters.getRmsdThreshold(), 2);
340
341                List<Point3d> centers = subunits.getOriginalCenters();
342                List<Integer> seqClusterId = subunits.getSequenceClusterIds();
343
344                List<Integer> permutations = new ArrayList<Integer>(centers.size());
345                double[] dSqs = new double[centers.size()];
346                boolean[] used = new boolean[centers.size()];
347                Arrays.fill(used, false);
348
349                for (int i = 0; i < centers.size(); i++) {
350                        Point3d tCenter = new Point3d(centers.get(i));
351                        transformation.transform(tCenter);
352                        int permutation = -1;
353                        double minDistSq = Double.MAX_VALUE;
354                        for (int j = 0; j < centers.size(); j++) {
355                                if (seqClusterId.get(i) == seqClusterId.get(j)) {
356                                        if (! used[j]) {
357                                                double dSq = tCenter.distanceSquared(centers.get(j));
358                                                if (dSq < minDistSq && dSq <= rmsdThresholdSq) {
359                                                        minDistSq = dSq;
360                                                        permutation = j;
361                                                        dSqs[j] = dSq;
362                                                }
363                                        }
364                                }
365                        }
366                        // can't map to itself
367                        if (permutations.size() == permutation) {
368                                permutation = -1;
369                        }
370
371                        if (permutation != -1) {
372                                used[permutation] = true;
373                        }
374
375
376                        permutations.add(permutation);
377                }
378
379                return permutations;
380        }
381
382
383
384        /**
385         * Returns the rise of a helix given the subunit centers of two adjacent
386         * subunits and the helix transformation
387         * @param transformation helix transformation
388         * @param p1 center of one subunit
389         * @param p2 center of an adjacent subunit
390         * @return
391         */
392        private static double getRise(Matrix4d transformation, Point3d p1, Point3d p2) {
393                AxisAngle4d axis = getAxisAngle(transformation);
394                Vector3d h = new Vector3d(axis.x, axis.y, axis.z);
395                Vector3d p = new Vector3d();
396                p.sub(p1, p2);
397                return p.dot(h);
398        }
399
400        /**
401         * Returns the pitch angle of the helix
402         * @param transformation helix transformation
403         * @return
404         */
405        private static double getAngle(Matrix4d transformation) {
406                return getAxisAngle(transformation).angle;
407        }
408
409        /**
410         * Returns the AxisAngle of the helix transformation
411         * @param transformation helix transformation
412         * @return
413         */
414        private static AxisAngle4d getAxisAngle(Matrix4d transformation) {
415                AxisAngle4d axis = new AxisAngle4d();
416                axis.set(transformation);
417                return axis;
418        }
419}