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/**
022 *
023 */
024package org.biojava.nbio.structure.symmetry.jmolScript;
025
026import org.biojava.nbio.structure.symmetry.axis.HelixAxisAligner;
027import org.biojava.nbio.structure.symmetry.core.QuatSymmetrySubunits;
028import org.jcolorbrewer.ColorBrewer;
029
030import javax.vecmath.*;
031
032import java.awt.*;
033import java.util.*;
034import java.util.List;
035
036
037/**
038 * @author Peter
039 *
040 */
041public class JmolSymmetryScriptGeneratorH extends JmolSymmetryScriptGenerator {
042        private static double AXIS_SCALE_FACTOR = 1.2;
043        private static double SIDE_CHAIN_EXTENSION = 6.0;
044        private HelixAxisAligner helixAxisAligner = null;
045        private String name = "";
046        private String defaultColoring = "";
047        private boolean onTheFly = false;
048
049        public JmolSymmetryScriptGeneratorH(HelixAxisAligner helixAxisAligner, String name) {
050                this.helixAxisAligner = helixAxisAligner;
051                this.name = name;
052        }
053
054        @Override
055        public void setOnTheFly(boolean onTheFly) {
056                this.onTheFly = onTheFly;
057        }
058
059        @Override
060        public int getZoom() {
061                int zoom = (int) Math.round(100.0/AXIS_SCALE_FACTOR);
062                return zoom;
063        }
064
065        /**
066         * Returns a Jmol script to set the default orientation for a structure
067         * @return Jmol script
068         */
069        @Override
070        public String getDefaultOrientation() {
071                StringBuilder s = new StringBuilder();
072                s.append(setCentroid());
073
074                Quat4d q = new Quat4d();
075                q.set(helixAxisAligner.getRotationMatrix());
076
077                // set orientation
078                s.append("moveto 0 quaternion{");
079                s.append(jMolFloat(q.x));
080                s.append(",");
081                s.append(jMolFloat(q.y));
082                s.append(",");
083                s.append(jMolFloat(q.z));
084                s.append(",");
085                s.append(jMolFloat(q.w));
086                s.append("};");
087                return s.toString();
088        }
089
090        @Override
091        public int getOrientationCount() {
092                return 4;
093        }
094
095        @Override
096        public String getOrientation(int index) {
097                StringBuilder s = new StringBuilder();
098                s.append(setCentroid());
099
100                // calculate  orientation
101                Matrix4d matrix = new Matrix4d();
102                matrix.setIdentity();
103                matrix.setRotation(new AxisAngle4d(1,0,0,Math.PI/2*(index+1)));
104                Quat4d r = new Quat4d();
105                r.set(new AxisAngle4d(1,0,0,index*Math.PI/2));
106                Quat4d q = new Quat4d();
107                q.set(helixAxisAligner.getRotationMatrix());
108                r.mul(q);
109
110                // set orientation
111                s.append("moveto 4 quaternion{");
112                s.append(jMolFloat(r.x));
113                s.append(",");
114                s.append(jMolFloat(r.y));
115                s.append(",");
116                s.append(jMolFloat(r.z));
117                s.append(",");
118                s.append(jMolFloat(r.w));
119                s.append("}");
120                s.append(";");
121                return s.toString();
122        }
123
124        /**
125         * Returns a Jmol script that sets a specific orientation and zoom
126         * to draw either axes or polyhedron
127         * @param index orientation index
128         * @return Jmol script
129         */
130        @Override
131        public String getOrientationWithZoom(int index) {
132                StringBuilder s = new StringBuilder();
133                s.append(getOrientation(index));
134                s.insert(s.length()-1, getZoom());
135                return s.toString();
136
137        }
138        /**
139         * Returns the name of a specific orientation
140         * @param index orientation index
141         * @return name of orientation
142         */
143        @Override
144        public String getOrientationName(int index) {
145                switch (index) {
146                        case 0: return "Side";
147                        case 1: return "Front";
148                        case 2: return "Other side";
149                        case 3: return "Back";
150                        default: return "";
151                }
152        }
153
154        /* (non-Javadoc)
155         * @see org.biojava.nbio.structure.quaternary.jmolScript.JMolSymmetryScriptInterface#getTransformation()
156         */
157        @Override
158        public Matrix4d getTransformation() {
159                return helixAxisAligner.getTransformation();
160        }
161
162
163        @Override
164        public void setDefaultColoring(String colorScript) {
165                this.defaultColoring = colorScript;
166        }
167
168        /**
169         * Returns a Jmol script that draws an invisible polyhedron around a structure.
170         * Use showPolyhedron() and hidePolyhedron() to toggle visibility.
171         * @return Jmol script
172         */
173        @Override
174        public String drawPolyhedron() {
175                StringBuilder s = new StringBuilder();
176
177                // for now, return empty script until this has been implemented properly
178//              if (s.length() == 0)
179//                      return s.toString();
180
181//              Point3d[] vertices = getRepeatUnitCenters();
182                List<Point3d> vertices = helixAxisAligner.getSubunits().getOriginalCenters();
183                vertices = extendUnitCenters(vertices);
184
185                int index = 0;
186                Color4f c = new Color4f(Color.MAGENTA);
187//              double width = getMaxExtension()*0.015;
188                double width = getMaxExtension()*0.007;
189
190//              List<List<Integer>> layerLines = helixAxisAligner.getHelixLayers().getByLargestContactsNotLowestAngle().getLayerLines();
191//              layerLines.addAll(helixAxisAligner.getHelixLayers().getByLowestAngle().getLayerLines());
192                List<List<Integer>> layerLines = helixAxisAligner.getHelixLayers().getByLargestContacts().getLayerLines();
193
194
195                for (List<Integer> line : layerLines) {
196                        s.append("draw polyhedron");
197                        s.append(name);
198                        s.append(index++);
199                        s.append(" line");
200                        for (int i: line) {
201                                s.append(getJmolPoint(vertices.get(i)));
202                        }
203                        s.append("width ");
204                        s.append(fDot2(width));
205                        s.append(" color");
206                        s.append(getJmolColor(c));
207                        s.append(" off;");
208                }
209
210                List<Point3d> interiorVertices = interiorCenters(vertices);
211
212                for (int i = 0; i < vertices.size(); i++) {
213                        s.append("draw polyhedron");
214                        s.append(name);
215                        s.append(index++);
216                        s.append(" line");
217                        s.append(getJmolPoint(vertices.get(i)));
218                        s.append(getJmolPoint(interiorVertices.get(i)));
219                        s.append("width ");
220                        s.append(fDot2(width));
221                        s.append(" color");
222                        s.append(getJmolColor(c));
223                        s.append(" off;");
224                }
225
226//              Matrix4d transformation = helixAxisTransformation.getHelixLayers().getByLowestAngle().getTransformation();
227//
228//              Point3d[] vertices2 = SuperPosition.clonePoint3dArray(vertices);
229//              for (Point3d p: vertices2) {
230//                      transformation.transform(p);
231//              }
232//
233//              // extend grid by one turn in +z direction
234//              for (int[] lineLoop: getLayerLines()) {
235//                      s.append("draw line");
236//                      s.append(name);
237//                      s.append(index++);
238//                      s.append(" line");
239//                      int count = 0;
240//                      for (int i: lineLoop) {
241//                              for (Point3d v: vertices) {
242//                                      if (v.distance(vertices2[i]) > 0.1) {
243//                                              count++;
244//                                              break;
245//                                      }
246//                              }
247//                              if (count > 0) {
248//                                      s.append(getJmolPoint(vertices2[i]));
249//                              } else {
250//                                      s.append(" ");
251//                              }
252//                      }
253//                      s.append("width ");
254//                  s.append(fDot2(width));
255//                      s.append(" color");
256//                      s.append(getJmolColor(c));
257//                      s.append(" off;");
258//              }
259//
260//
261//              // extend grid by one turn into -z direction
262//              Point3d[] vertices3 = SuperPosition.clonePoint3dArray(vertices);
263//              Matrix4d m = new Matrix4d();
264//              m.set(transformation);
265//              m.invert();
266//              for (Point3d p: vertices3) {
267//                      m.transform(p);
268//              }
269//
270//              for (int[] lineLoop: getLayerLines()) {
271//                      s.append("draw line");
272//                      s.append(name);
273//                      s.append(index++);
274//                      s.append(" line");
275//                      int count = 0;
276//                      for (int i: lineLoop) {
277//                              for (Point3d v: vertices) {
278//                                      if (v.distance(vertices3[i]) > 0.1) {
279//                                              count++;
280//                                      }
281//                              }
282//                              if (count > 0) {
283//                                      s.append(getJmolPoint(vertices3[i]));
284//                              } else {
285//                                      s.append(" ");
286//                              }
287//                      }
288//                      s.append("width ");
289//                  s.append(fDot2(width));
290//                      s.append(" color");
291//                      s.append(getJmolColor(c));
292//                      s.append(" off;");
293//              }
294
295                return s.toString();
296        }
297
298        @Override
299        public String hidePolyhedron() {
300                return "draw polyhedron* off;";
301        }
302
303        @Override
304        public String showPolyhedron() {
305                return "draw polyhedron* on;";
306        }
307
308        /**
309         * Returns a Jmol script that draws symmetry or inertia axes for a structure.
310         * Use showAxes() and hideAxes() to toggle visibility.
311         * @return Jmol script
312         */
313        @Override
314        public String drawAxes() {
315                StringBuilder s = new StringBuilder();
316//              Point3d centroid = helixAxisAligner.getCentroid();
317//              System.out.println("Centroid: " + helixAxisAligner.getCentroid());
318//              Point3d centroid = helixAxisAligner.getGeometricCenter();
319                Point3d centroid = helixAxisAligner.calcCenterOfRotation();
320//              System.out.println("Geometric center: " + centroid);
321                AxisAngle4d axisAngle = helixAxisAligner.getHelixLayers().getByLowestAngle().getAxisAngle();
322                Vector3d axis = new Vector3d(axisAngle.x, axisAngle.y, axisAngle.z);
323
324                s.append("draw axesHelical");
325                s.append(name);
326                s.append(0);
327                s.append(" ");
328                s.append("line");
329                Point3d v1 = new Point3d(axis);
330                v1.scale(AXIS_SCALE_FACTOR*(helixAxisAligner.getDimension().y + SIDE_CHAIN_EXTENSION));
331                Point3d v2 = new Point3d(v1);
332                v2.negate();
333                v1.add(centroid);
334                v2.add(centroid);
335                s.append(getJmolPoint(v1));
336                s.append(getJmolPoint(v2));
337                s.append("width 1.0 ");
338                s.append(" color red");
339                s.append(" off;");
340                return s.toString();
341        }
342
343        /**
344         * Returns a Jmol script to hide axes
345         * @return Jmol script
346         */
347        @Override
348        public String hideAxes() {
349                return "draw axes* off;";
350        }
351
352        /**
353         * Returns a Jmol script to show axes
354         * @return Jmol script
355         */
356        @Override
357        public String showAxes() {
358                return "draw axes* on;";
359        }
360
361        /**
362         * Returns a Jmol script that displays a symmetry polyhedron and symmetry axes
363         * and then loop through different orientations
364         * @return Jmol script
365         */
366        @Override
367        public String playOrientations() {
368                StringBuilder s = new StringBuilder();
369
370                // draw footer
371                s.append(drawFooter("Symmetry Helical", "white"));
372
373                // draw polygon
374                s.append(drawPolyhedron()); // draw invisibly
375                s.append(showPolyhedron());
376
377                // draw axes
378                s.append(drawAxes());
379                s.append(showAxes());
380
381                // loop over all orientations with 4 sec. delay
382                for (int i = 0; i < getOrientationCount(); i++) {
383                        s.append(deleteHeader());
384                        s.append(getOrientationWithZoom(i));
385                        s.append(drawHeader(getOrientationName(i), "white"));
386                        s.append("delay 4;");
387                }
388
389                // go back to first orientation
390                s.append(deleteHeader());
391                s.append(getOrientationWithZoom(0));
392                s.append(drawHeader(getOrientationName(0), "white"));
393
394                return s.toString();
395        }
396
397
398        /**
399         * Returns a Jmol script that colors the subunits of a structure by different colors
400         * @return
401         */
402        @Override
403        public String colorBySubunit() {
404                QuatSymmetrySubunits subunits = helixAxisAligner.getSubunits();
405                List<Integer> modelNumbers = subunits.getModelNumbers();
406                List<String> chainIds = subunits.getChainIds();
407                List<List<Integer>> orbits = helixAxisAligner.getOrbits();
408
409                Color[] col = ColorBrewer.Spectral.getColorPalette(orbits.size());
410                Color4f[] colors = ColorConverter.convertColor4f(col);
411                int half = colors.length/2;
412                for (int i = 0; i < half; i++) {
413                        if (i % 2 != 0) {
414                           Color4f temp = colors[i];
415                           colors[i] = colors[half+i];
416                           colors[half+i] = temp;
417                        }
418                }
419                Map<Color4f, List<String>> colorMap = new HashMap<Color4f, List<String>>();
420
421                for (int i = 0; i < orbits.size(); i++) {
422                                for (Integer j: orbits.get(i)) {
423                                        Color4f c = colors[i];
424                                        List<String> ids = colorMap.get(c);
425                                        if (ids == null) {
426                                                ids = new ArrayList<String>();
427                                                colorMap.put(c,  ids);
428                                        }
429                                        String id = getChainSpecification(modelNumbers, chainIds, j);
430                                        ids.add(id);
431                                }
432                }
433                String coloring = defaultColoring + getJmolColorScript(colorMap);
434                return coloring;
435        }
436
437        /**
438         * Returns a Jmol script that colors subunits by their sequence cluster ids.
439         * @return Jmol script
440         */
441        @Override
442        public String colorBySequenceCluster() {
443                QuatSymmetrySubunits subunits = helixAxisAligner.getSubunits();
444                int n = subunits.getSubunitCount();
445                List<Integer> modelNumbers = subunits.getModelNumbers();
446                List<String> chainIds = subunits.getChainIds();
447                List<Integer> seqClusterIds = subunits.getClusterIds();
448                int clusters = Collections.max(seqClusterIds) + 1;
449                Color[] col = ColorBrewer.BrBG.getColorPalette(clusters);
450                Color4f[] colors = ColorConverter.convertColor4f(col);
451                Map<Color4f, List<String>> colorMap = new HashMap<Color4f, List<String>>();
452
453                for (int i = 0; i < n; i++) {
454                        Color4f c = colors[seqClusterIds.get(i)];
455                        List<String> ids = colorMap.get(c);
456                        if (ids == null) {
457                                ids = new ArrayList<String>();
458                                colorMap.put(c,  ids);
459                        }
460                        String id = getChainSpecification(modelNumbers, chainIds, i);
461                        ids.add(id);
462
463                }
464                String coloring = defaultColoring + getJmolColorScript(colorMap);
465                return coloring;
466        }
467
468
469        /**
470         * Returns a Jmol script that colors subunits to highlight the symmetry within a structure
471         * Different subunits should have a consistent color scheme or different shade of the same colors
472         * @return Jmol script
473         */
474        @Override
475        public String colorBySymmetry() {
476                List<List<Integer>> units = helixAxisAligner.getHelixLayers().getByLargestContacts().getLayerLines();
477                units = orientLayerLines(units);
478                QuatSymmetrySubunits subunits = helixAxisAligner.getSubunits();
479                List<Integer> modelNumbers = subunits.getModelNumbers();
480                List<String> chainIds = subunits.getChainIds();
481                List<Integer> clusterIds = subunits.getClusterIds();
482                int clusterCount = Collections.max(clusterIds) + 1;
483
484                Map<Color4f, List<String>> colorMap = new HashMap<Color4f, List<String>>();
485
486                int maxLen = 0;
487                for (List<Integer> unit: units) {
488                        maxLen = Math.max(maxLen,  unit.size());
489                }
490
491//              Color4f[] colors = getSymmetryColors(permutation.size());
492                Color4f[] colors = getSymmetryColors(subunits.getSubunitCount());
493                int count = 0;
494
495                for (int i = 0; i < maxLen; i++) {
496                        for (int j = 0; j < units.size(); j++) {
497                                int m = units.get(j).size();
498                                if (i < m) {
499                                        int subunit = units.get(j).get(i);
500                                        int cluster = clusterIds.get(subunit);
501                                        float scale = 0.3f + 0.7f * (cluster+1)/clusterCount;
502                                        Color4f c = new Color4f(colors[count]);
503                                        count++;
504                                        c.scale(scale);
505                                        List<String> ids = colorMap.get(c);
506                                        if (ids == null) {
507                                                ids = new ArrayList<String>();
508                                                colorMap.put(c,  ids);
509                                        }
510                                        String id = getChainSpecification(modelNumbers, chainIds, subunit);
511                                        ids.add(id);
512                                }
513                        }
514                }
515
516                String coloring = defaultColoring + getJmolColorScript(colorMap);
517                return coloring;
518        }
519
520        private String getChainSpecification(List<Integer> modelNumbers, List<String> chainIds, int subunit) {
521                if (onTheFly) {
522                        if (Collections.max(modelNumbers) > 1) {
523                                return chainIds.get(subunit) + "&symop=" + (modelNumbers.get(subunit)+1);
524                        } else {
525                                // if there is only a single symop, Jmol does not accept the symop syntax
526                                return chainIds.get(subunit);
527                        }
528                } else {
529                        return chainIds.get(subunit) + "/" + (modelNumbers.get(subunit)+1);
530                }
531        }
532
533        /**
534         * Orients layer lines from lowest y-axis value to largest y-axis value
535         */
536        private List<List<Integer>> orientLayerLines(List<List<Integer>> layerLines) {
537                Matrix4d transformation = helixAxisAligner.getTransformation();
538                List<Point3d> centers = helixAxisAligner.getSubunits().getOriginalCenters();
539
540                for (int i = 0; i < layerLines.size(); i++) {
541                        List<Integer> layerLine = layerLines.get(i);
542
543                        // get center of first subunit in layerline and transform to standard orientation (helix axis aligned with y-axis)
544                        int first = layerLine.get(0);
545                        Point3d firstSubunit = new Point3d(centers.get(first));
546                        transformation.transform(firstSubunit);
547
548                        // get center of last subunit in layerline and transform to standard orientation (helix axis aligned with y-axis)
549                        int last = layerLine.get(layerLine.size()-1);
550                        Point3d lastSubunit = new Point3d(centers.get(last));
551                        transformation.transform(lastSubunit);
552
553                        // a layerline should start at the lowest y-value, so all layerlines have a consistent direction from -y value to +y value
554                        if (firstSubunit.y > lastSubunit.y) {
555//                              System.out.println("reorienting layer line: " + layerLine);
556                                Collections.reverse(layerLine);
557                        }
558                }
559                return layerLines;
560        }
561
562        // --- protected methods ---
563        /**
564         * Returns the maximum extension (length) of structure
565         * @return
566         */
567        protected double getMaxExtension() {
568                Vector3d dimension = helixAxisAligner.getDimension();
569                double maxExtension = Math.max(dimension.x, dimension.y);
570                maxExtension = Math.max(maxExtension, dimension.z);
571                return maxExtension;
572        }
573
574        private List<Point3d> extendUnitCenters(List<Point3d> points) {
575                Matrix4d transformation = helixAxisAligner.getTransformation();
576                Matrix4d reversetransformation = helixAxisAligner.getReverseTransformation();
577                double radius = helixAxisAligner.getRadius();
578                radius += SIDE_CHAIN_EXTENSION; // add space for side chains
579//              System.out.println("XZRadius: " + radius);
580
581
582                List<Point3d> ePoints = new ArrayList<Point3d>(points.size());
583                for (Point3d p: points) {
584                        Point3d q = new Point3d(p);
585                        transformation.transform(q);
586                        double d = Math.sqrt(q.x*q.x + q.z*q.z);
587                        double scale = radius/d;
588                        q.x *= scale;
589                        q.z *= scale;
590                        reversetransformation.transform(q);
591                        ePoints.add(q);
592                }
593                return ePoints;
594        }
595
596        private List<Point3d> interiorCenters(List<Point3d> points) {
597                Matrix4d transformation = helixAxisAligner.getTransformation();
598                Matrix4d reversetransformation = helixAxisAligner.getReverseTransformation();
599
600                List<Point3d> ePoints = new ArrayList<Point3d>(points.size());
601                for (Point3d p: points) {
602                        Point3d q = new Point3d(p);
603                        transformation.transform(q);
604                        q.x = 0;
605                        q.z = 0;
606                        reversetransformation.transform(q);
607                        ePoints.add(q);
608                }
609                return ePoints;
610        }
611
612        private String drawHeader(String text, String color) {
613                StringBuilder s = new StringBuilder();
614                s.append("set echo top center;");
615                s.append("color echo ");
616                s.append(color);
617                s.append(";");
618                s.append("font echo 24 sanserif;");
619                s.append("echo ");
620                s.append(text);
621                s.append(";");
622                return s.toString();
623        }
624
625        private String deleteHeader() {
626                return "set echo top center;echo ;";
627        }
628
629        private String drawFooter(String text, String color) {
630                StringBuilder s = new StringBuilder();
631                s.append("set echo bottom center;");
632                s.append("color echo ");
633                s.append(color);
634                s.append(";");
635                s.append("font echo 24 sanserif;");
636                s.append("echo ").append(text);
637                s.append(";");
638                return s.toString();
639        }
640
641        /**
642         * Returns a unique color palette based on point group
643         * @param nColors
644         * @return
645         */
646        private Color4f[] getSymmetryColors(int nColors) {
647                int offset = 0;
648                int dMax = nColors + offset;
649                Color[] col = ColorBrewer.Spectral.getColorPalette(dMax);
650                Color4f[] colors = ColorConverter.convertColor4f(col);
651                System.arraycopy(colors, offset, colors, 0, dMax-offset);
652                return colors;
653        }
654
655        private String setCentroid() {
656                // calculate center of rotation
657        //      Point3d centroid = axisTransformation.getGeometricCenter();
658        //      Point3d centroid = helixAxisAligner.getCentroid();
659                Point3d centroid = helixAxisAligner.calcCenterOfRotation();
660
661                // set centroid
662                StringBuilder s = new StringBuilder();
663                s.append("center");
664                s.append(getJmolPoint(centroid));
665                s.append(";");
666                return s.toString();
667        }
668
669}