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 * Created on May 21, 2006
021 *
022 */
023package org.biojava.nbio.structure.align;
024
025import org.biojava.nbio.structure.*;
026import org.biojava.nbio.structure.align.ce.GuiWrapper;
027import org.biojava.nbio.structure.align.helper.AlignUtils;
028import org.biojava.nbio.structure.align.helper.JointFragments;
029import org.biojava.nbio.structure.align.pairwise.*;
030import org.biojava.nbio.structure.geometry.Matrices;
031import org.biojava.nbio.structure.geometry.SuperPositions;
032import org.biojava.nbio.structure.io.PDBFileParser;
033import org.biojava.nbio.structure.io.PDBFileReader;
034import org.biojava.nbio.structure.jama.Matrix;
035import org.slf4j.Logger;
036import org.slf4j.LoggerFactory;
037
038import java.io.FileOutputStream;
039import java.io.InputStream;
040import java.io.PrintStream;
041import java.util.ArrayList;
042import java.util.Collections;
043import java.util.Comparator;
044import java.util.List;
045
046import javax.vecmath.Matrix4d;
047
048/**
049 * Perform a pairwise protein structure superimposition.
050 *
051 * <p>
052 * The algorithm is a distance matrix based, rigid body protein structure
053 * superimposition. It is based on a variation of the PSC++ algorithm provided
054 * by Peter Lackner (Peter.Lackner@sbg.ac.at, personal communication) .
055 * </p>
056 *
057 *
058 *
059 * <h2>Example</h2>
060 *
061 * <pre>
062 *  public void run(){
063 *
064 *              // first load two example structures
065 *              {@link InputStream} inStream1 = this.getClass().getResourceAsStream("/files/5pti.pdb");
066 *              {@link InputStream} inStream2 = this.getClass().getResourceAsStream("/files/1tap.pdb");
067 *
068 *              {@link Structure} structure1 = null;
069 *              {@link Structure} structure2 = null;
070 *
071 *              {@link PDBFileParser} pdbpars = new {@link PDBFileParser}();
072 *              structure1 = pdbpars.parsePDBFile(inStream1) ;
073 *              structure2 = pdbpars.parsePDBFile(inStream2);
074 *
075 *
076 *              // calculate structure superimposition for two complete structures
077 *              {@link StructurePairAligner} aligner = new {@link StructurePairAligner}();
078 *
079 *
080 *                      // align the full 2 structures with default parameters.
081 *                      // see StructurePairAligner for more options and how to align
082 *                      // any set of Atoms
083 *                      aligner.align(structure1,structure2);
084 *
085 *                      {@link AlternativeAlignment}[] aligs = aligner.getAlignments();
086 *                      {@link AlternativeAlignment} a = aligs[0];
087 *                      System.out.println(a);
088 *
089 *                      //display the alignment in Jmol
090 *
091 *                      // first get an artificial structure for the alignment
092 *                      {@link Structure} artificial = a.getAlignedStructure(structure1, structure2);
093 *
094 *
095 *                      // and then send it to Jmol (only will work if Jmol is in the Classpath)
096 *
097 *                      BiojavaJmol jmol = new BiojavaJmol();
098 *                      jmol.setTitle(artificial.getName());
099 *                      jmol.setStructure(artificial);
100 *
101 *                      // color the two structures
102 *
103 *
104 *                      jmol.evalString("select *; backbone 0.4; wireframe off; spacefill off; " +
105 *                                      "select not protein and not solvent; spacefill on;");
106 *                      jmol.evalString("select *"+"/1 ; color red; model 1; ");
107 *
108 *
109 *                      // now color the equivalent residues ...
110 *
111 *                      String[] pdb1 = a.getPDBresnum1();
112 *                      for (String res : pdb1 ){
113 *                              jmol.evalString("select " + res + "/1 ; backbone 0.6; color white;");
114 *                      }
115 *
116 *                      jmol.evalString("select *"+"/2; color blue; model 2;");
117 *                      String[] pdb2 = a.getPDBresnum2();
118 *                      for (String res :pdb2 ){
119 *                              jmol.evalString("select " + res + "/2 ; backbone 0.6; color yellow;");
120 *                      }
121 *
122 *
123 *                      // now show both models again.
124 *                      jmol.evalString("model 0;");
125 *
126 *      }
127 * </pre>
128 *
129 *
130 *
131 * @author Andreas Prlic
132 * @author Peter Lackner
133 * @since 1.4
134 * @version %I% %G%
135 */
136public class StructurePairAligner {
137
138        private final static Logger logger = LoggerFactory
139                        .getLogger(StructurePairAligner.class);
140
141        private AlternativeAlignment[] alts;
142        private Matrix distanceMatrix;
143        private StrucAligParameters params;
144        private FragmentPair[] fragPairs;
145
146        private final List<AlignmentProgressListener> listeners = new ArrayList<>();
147
148        public StructurePairAligner() {
149                super();
150                params = StrucAligParameters.getDefaultParameters();
151                reset();
152                alts = new AlternativeAlignment[0];
153                distanceMatrix = new Matrix(0, 0);
154        }
155
156        public void addProgressListener(AlignmentProgressListener li) {
157                listeners.add(li);
158        }
159
160        public void clearListeners() {
161                listeners.clear();
162        }
163
164        /**
165         * example usage of this class
166         *
167         * @param args
168         */
169        public static void main(String[] args) throws Exception {
170                // UPDATE THE FOLLOWING LINES TO MATCH YOUR SETUP
171
172                PDBFileReader pdbr = new PDBFileReader();
173                pdbr.setPath("/Users/andreas/WORK/PDB/");
174
175                // String pdb1 = "1crl";
176                // String pdb2 = "1ede";
177
178                String pdb1 = "1buz";
179                String pdb2 = "1ali";
180                String outputfile = "/tmp/alig_" + pdb1 + "_" + pdb2 + ".pdb";
181
182                // NO NEED TO DO CHANGE ANYTHING BELOW HERE...
183
184                StructurePairAligner sc = new StructurePairAligner();
185
186                // step1 : read molecules
187
188                logger.info("aligning {} vs. {}", pdb1, pdb2);
189
190                Structure s1 = pdbr.getStructureById(pdb1);
191                Structure s2 = pdbr.getStructureById(pdb2);
192
193                // step 2 : do the calculations
194                sc.align(s1, s2);
195
196                AlternativeAlignment[] aligs = sc.getAlignments();
197
198                // cluster similar results together
199                ClusterAltAligs.cluster(aligs);
200
201                // print the result:
202                // the AlternativeAlignment object gives also access to rotation
203                // matrices / shift vectors.
204                for (AlternativeAlignment aa : aligs) {
205                        logger.info("Alternative Alignment: {}", aa);
206                }
207
208                // convert AlternativeAlignemnt 1 to PDB file, so it can be opened with
209                // a viewer (e.g. Jmol, Rasmol)
210
211                if (aligs.length > 0) {
212                        AlternativeAlignment aa1 = aligs[0];
213                        String pdbstr = aa1.toPDB(s1, s2);
214
215                        logger.info("writing alignment to {}", outputfile);
216                        FileOutputStream out = new FileOutputStream(outputfile);
217                        PrintStream p = new PrintStream(out);
218
219                        p.println(pdbstr);
220
221                        p.close();
222                        out.close();
223                }
224
225                // display the alignment in Jmol
226                // only will work if Jmol is in the Classpath
227
228                if (aligs.length > 0) {
229
230                        if (!GuiWrapper.isGuiModuleInstalled()) {
231                                logger.error("Could not find structure-gui modules in classpath, please install first!");
232                        }
233
234                }
235
236        }
237
238        private void reset() {
239                alts = new AlternativeAlignment[0];
240                distanceMatrix = new Matrix(0, 0);
241                fragPairs = new FragmentPair[0];
242
243        }
244
245        /**
246         * get the results of step 1 - the FragmentPairs used for seeding the
247         * alignment
248         *
249         * @return a FragmentPair[] array
250         */
251
252        public FragmentPair[] getFragmentPairs() {
253                return fragPairs;
254        }
255
256        public void setFragmentPairs(FragmentPair[] fragPairs) {
257                this.fragPairs = fragPairs;
258        }
259
260        /**
261         * return the alternative alignments that can be found for the two
262         * structures
263         *
264         * @return AlternativeAlignment[] array
265         */
266        public AlternativeAlignment[] getAlignments() {
267                return alts;
268        }
269
270        /**
271         * return the difference of distance matrix between the two structures
272         *
273         * @return a Matrix
274         */
275        public Matrix getDistMat() {
276                return distanceMatrix;
277        }
278
279        /**
280         * get the parameters.
281         *
282         * @return the Parameters.
283         */
284        public StrucAligParameters getParams() {
285                return params;
286        }
287
288        /**
289         * set the parameters to be used for the algorithm
290         *
291         * @param params
292         *            the Parameter object
293         */
294        public void setParams(StrucAligParameters params) {
295                this.params = params;
296        }
297
298        /**
299         * Calculate the alignment between the two full structures with default
300         * parameters
301         *
302         * @param s1
303         * @param s2
304         * @throws StructureException
305         */
306        public void align(Structure s1, Structure s2) throws StructureException {
307
308                align(s1, s2, params);
309        }
310
311        /**
312         * Calculate the alignment between the two full structures with user
313         * provided parameters
314         *
315         * @param s1
316         * @param s2
317         * @param params
318         * @throws StructureException
319         */
320        public void align(Structure s1, Structure s2, StrucAligParameters params)
321                        throws StructureException {
322                // step 1 convert the structures to Atom Arrays
323
324                Atom[] ca1 = getAlignmentAtoms(s1);
325                Atom[] ca2 = getAlignmentAtoms(s2);
326
327                notifyStartingAlignment(s1.getName(), ca1, s2.getName(), ca2);
328                align(ca1, ca2, params);
329        }
330
331        /**
332         * Align two chains from the structures. Uses default parameters.
333         *
334         * @param s1
335         * @param chainId1
336         * @param s2
337         * @param chainId2
338         */
339        public void align(Structure s1, String chainId1, Structure s2,
340                        String chainId2) throws StructureException {
341                align(s1, chainId1, s2, chainId2, params);
342        }
343
344        /**
345         * Aligns two chains from the structures using user provided parameters.
346         *
347         * @param s1
348         * @param chainId1
349         * @param s2
350         * @param chainId2
351         * @param params
352         * @throws StructureException
353         */
354        public void align(Structure s1, String chainId1, Structure s2,
355                        String chainId2, StrucAligParameters params)
356                        throws StructureException {
357                reset();
358                this.params = params;
359
360                Chain c1 = s1.getPolyChainByPDB(chainId1);
361                Chain c2 = s2.getPolyChainByPDB(chainId2);
362
363                Structure s3 = new StructureImpl();
364                s3.addChain(c1);
365
366                Structure s4 = new StructureImpl();
367                s4.addChain(c2);
368
369                Atom[] ca1 = getAlignmentAtoms(s3);
370                Atom[] ca2 = getAlignmentAtoms(s4);
371
372                notifyStartingAlignment(s1.getName(), ca1, s2.getName(), ca2);
373                align(ca1, ca2, params);
374        }
375
376        /**
377         * Returns the atoms that are being used for the alignment. (E.g. Calpha
378         * only, etc.)
379         *
380         * @param s
381         * @return an array of Atoms objects
382         */
383        public Atom[] getAlignmentAtoms(Structure s) {
384                String[] atomNames = params.getUsedAtomNames();
385                return StructureTools.getAtomArray(s, atomNames);
386        }
387
388        /**
389         * Calculate the protein structure superimposition, between two sets of
390         * atoms.
391         *
392         * @param ca1
393         *            set of Atoms of structure 1
394         * @param ca2
395         *            set of Atoms of structure 2
396         * @param params
397         *            the parameters to use for the alignment
398         * @throws StructureException
399         */
400        public void align(Atom[] ca1, Atom[] ca2, StrucAligParameters params)
401                        throws StructureException {
402
403                reset();
404                this.params = params;
405
406                long timeStart = System.currentTimeMillis();
407
408                // step 1 get all Diagonals of length X that are similar between both
409                // structures
410                logger.debug(" length atoms1:" + ca1.length);
411                logger.debug(" length atoms2:" + ca2.length);
412
413                logger.debug("step 1 - get fragments with similar intramolecular distances ");
414
415                int k = params.getDiagonalDistance();
416                int k2 = params.getDiagonalDistance2();
417                int fragmentLength = params.getFragmentLength();
418
419                if (ca1.length < (fragmentLength + 1)) {
420                        throw new StructureException("structure 1 too short (" + ca1.length
421                                        + "), can not align");
422                }
423                if (ca2.length < (fragmentLength + 1)) {
424                        throw new StructureException("structure 2 too short (" + ca2.length
425                                        + "), can not align");
426                }
427                int rows = ca1.length - fragmentLength + 1;
428                int cols = ca2.length - fragmentLength + 1;
429                distanceMatrix = new Matrix(rows, cols, 0.0);
430
431                double[] dist1 = AlignUtils.getDiagonalAtK(ca1, k);
432
433                double[] dist2 = AlignUtils.getDiagonalAtK(ca2, k);
434                double[] dist3 = new double[0];
435                double[] dist4 = new double[0];
436                if (k2 > 0) {
437                        dist3 = AlignUtils.getDiagonalAtK(ca1, k2);
438                        dist4 = AlignUtils.getDiagonalAtK(ca2, k2);
439                }
440
441                double[][] utmp = new double[][] { { 0, 0, 1 } };
442                Atom unitvector = new AtomImpl();
443                unitvector.setCoords(utmp[0]);
444
445                List<FragmentPair> fragments = new ArrayList<>();
446
447                for (int i = 0; i < rows; i++) {
448
449                        Atom[] catmp1 = AlignUtils.getFragment(ca1, i, fragmentLength);
450                        Atom center1 = AlignUtils.getCenter(ca1, i, fragmentLength);
451
452                        for (int j = 0; j < cols; j++) {
453
454                                double rdd1 = AlignUtils.rms_dk_diag(dist1, dist2, i, j,
455                                                fragmentLength, k);
456                                double rdd2 = 0;
457                                if (k2 > 0)
458                                        rdd2 = AlignUtils.rms_dk_diag(dist3, dist4, i, j,
459                                                        fragmentLength, k2);
460                                double rdd = rdd1 + rdd2;
461                                distanceMatrix.set(i, j, rdd);
462
463                                if (rdd < params.getFragmentMiniDistance()) {
464                                        FragmentPair f = new FragmentPair(fragmentLength, i, j);
465                                        Atom[] catmp2 = AlignUtils.getFragment(ca2, j,
466                                                        fragmentLength);
467                                        Atom center2 = AlignUtils.getCenter(ca2, j,
468                                                        fragmentLength);
469
470                                        f.setCenter1(center1);
471                                        f.setCenter2(center2);
472
473                                        Matrix4d t = SuperPositions.superpose(
474                                                        Calc.atomsToPoints(catmp1),
475                                                        Calc.atomsToPoints(catmp2));
476
477                                        Matrix rotmat = Matrices.getRotationJAMA(t);
478                                        f.setRot(rotmat);
479
480                                        Atom aunitv = (Atom) unitvector.clone();
481                                        Calc.rotate(aunitv, rotmat);
482                                        f.setUnitv(aunitv);
483
484                                        boolean doNotAdd = false;
485                                        if (params.reduceInitialFragments()) {
486                                                doNotAdd = FragmentJoiner.reduceFragments(
487                                                                fragments, f, distanceMatrix);
488
489                                        }
490                                        if (doNotAdd)
491                                                continue;
492
493                                        fragments.add(f);
494                                }
495                        }
496                }
497
498                notifyFragmentListeners(fragments);
499
500                FragmentPair[] fp = fragments.toArray(new FragmentPair[0]);
501                setFragmentPairs(fp);
502
503                logger.debug(" got # fragment pairs: {}", fp.length);
504
505                logger.debug("step 2 - join fragments");
506
507                // step 2 combine them to possible models
508                FragmentJoiner joiner = new FragmentJoiner();
509
510                JointFragments[] frags;
511
512                if (params.isJoinFast()) {
513                        // apply the quick alignment procedure.
514                        // less quality in alignments, better for DB searches...
515                        frags = joiner.approach_ap3(ca1, ca2, fp, params);
516
517                        joiner.extendFragments(ca1, ca2, frags, params);
518
519                } else if (params.isJoinPlo()) {
520                        // this approach by StrComPy (peter lackner):
521                        frags = joiner.frag_pairwise_compat(fp, params.getAngleDiff(),
522                                        params.getFragCompat(), params.getMaxrefine());
523
524                } else {
525
526                        // my first implementation
527                        frags = joiner.approach_ap3(ca1, ca2, fp, params);
528                }
529
530                notifyJointFragments(frags);
531
532                logger.debug(" number joint fragments: {}", frags.length);
533
534                logger.debug("step 3 - refine alignments");
535
536                List<AlternativeAlignment> aas = new ArrayList<>();
537                for (int i = 0; i < frags.length; i++) {
538                        JointFragments f = frags[i];
539                        AlternativeAlignment a = new AlternativeAlignment();
540                        a.apairs_from_idxlst(f);
541                        a.setAltAligNumber(i + 1);
542                        a.setDistanceMatrix(distanceMatrix);
543
544                        try {
545                                if (params.getMaxIter() > 0) {
546
547                                        a.refine(params, ca1, ca2);
548                                } else {
549
550                                        a.finish(params, ca1, ca2);
551
552                                }
553                        } catch (StructureException e) {
554                                logger.error("Refinement of fragment {} failed", i, e);
555                        }
556
557                        a.calcScores(ca1, ca2);
558                        aas.add(a);
559                }
560
561                // sort the alternative alignments
562                Comparator<AlternativeAlignment> comp = new AltAligComparator();
563                aas.sort(comp);
564                Collections.reverse(aas);
565
566                alts = aas.toArray(new AlternativeAlignment[0]);
567                // do final numbering of alternative solutions
568                int aanbr = 0;
569                for (AlternativeAlignment a : alts) {
570                        aanbr++;
571                        a.setAltAligNumber(aanbr);
572                }
573
574                logger.debug("total calculation time: {} ms.",
575                                (System.currentTimeMillis() - timeStart));
576        }
577
578        private void notifyStartingAlignment(String name1, Atom[] ca1,
579                        String name2, Atom[] ca2) {
580                for (AlignmentProgressListener li : listeners) {
581                        li.startingAlignment(name1, ca1, name2, ca2);
582                }
583        }
584
585        private void notifyFragmentListeners(List<FragmentPair> fragments) {
586
587                for (AlignmentProgressListener li : listeners) {
588                        li.calculatedFragmentPairs(fragments);
589                }
590
591        }
592
593        private void notifyJointFragments(JointFragments[] fragments) {
594                for (AlignmentProgressListener li : listeners) {
595                        li.jointFragments(fragments);
596                }
597        }
598
599}