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        AlternativeAlignment[] alts;
142        Matrix distanceMatrix;
143        StrucAligParameters params;
144        FragmentPair[] fragPairs;
145
146        List<AlignmentProgressListener> listeners = new ArrayList<AlignmentProgressListener>();
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         *
393         *
394         * @param ca1
395         *            set of Atoms of structure 1
396         * @param ca2
397         *            set of Atoms of structure 2
398         * @param params
399         *            the parameters to use for the alignment
400         * @throws StructureException
401         */
402        public void align(Atom[] ca1, Atom[] ca2, StrucAligParameters params)
403                        throws StructureException {
404
405                reset();
406                this.params = params;
407
408                long timeStart = System.currentTimeMillis();
409
410                // step 1 get all Diagonals of length X that are similar between both
411                // structures
412                logger.debug(" length atoms1:" + ca1.length);
413                logger.debug(" length atoms2:" + ca2.length);
414
415                logger.debug("step 1 - get fragments with similar intramolecular distances ");
416
417                int k = params.getDiagonalDistance();
418                int k2 = params.getDiagonalDistance2();
419                int fragmentLength = params.getFragmentLength();
420
421                if (ca1.length < (fragmentLength + 1)) {
422                        throw new StructureException("structure 1 too short (" + ca1.length
423                                        + "), can not align");
424                }
425                if (ca2.length < (fragmentLength + 1)) {
426                        throw new StructureException("structure 2 too short (" + ca2.length
427                                        + "), can not align");
428                }
429                int rows = ca1.length - fragmentLength + 1;
430                int cols = ca2.length - fragmentLength + 1;
431                distanceMatrix = new Matrix(rows, cols, 0.0);
432
433                double[] dist1 = AlignUtils.getDiagonalAtK(ca1, k);
434
435                double[] dist2 = AlignUtils.getDiagonalAtK(ca2, k);
436                double[] dist3 = new double[0];
437                double[] dist4 = new double[0];
438                if (k2 > 0) {
439                        dist3 = AlignUtils.getDiagonalAtK(ca1, k2);
440                        dist4 = AlignUtils.getDiagonalAtK(ca2, k2);
441                }
442
443                double[][] utmp = new double[][] { { 0, 0, 1 } };
444                Atom unitvector = new AtomImpl();
445                unitvector.setCoords(utmp[0]);
446
447                List<FragmentPair> fragments = new ArrayList<FragmentPair>();
448
449                for (int i = 0; i < rows; i++) {
450
451                        Atom[] catmp1 = AlignUtils.getFragment(ca1, i, fragmentLength);
452                        Atom center1 = AlignUtils.getCenter(ca1, i, fragmentLength);
453
454                        for (int j = 0; j < cols; j++) {
455
456                                double rdd1 = AlignUtils.rms_dk_diag(dist1, dist2, i, j,
457                                                fragmentLength, k);
458                                double rdd2 = 0;
459                                if (k2 > 0)
460                                        rdd2 = AlignUtils.rms_dk_diag(dist3, dist4, i, j,
461                                                        fragmentLength, k2);
462                                double rdd = rdd1 + rdd2;
463                                distanceMatrix.set(i, j, rdd);
464
465                                if (rdd < params.getFragmentMiniDistance()) {
466                                        FragmentPair f = new FragmentPair(fragmentLength, i, j);
467                                        Atom[] catmp2 = AlignUtils.getFragment(ca2, j,
468                                                        fragmentLength);
469                                        Atom center2 = AlignUtils.getCenter(ca2, j,
470                                                        fragmentLength);
471
472                                        f.setCenter1(center1);
473                                        f.setCenter2(center2);
474
475                                        Matrix4d t = SuperPositions.superpose(
476                                                        Calc.atomsToPoints(catmp1),
477                                                        Calc.atomsToPoints(catmp2));
478
479                                        Matrix rotmat = Matrices.getRotationJAMA(t);
480                                        f.setRot(rotmat);
481
482                                        Atom aunitv = (Atom) unitvector.clone();
483                                        Calc.rotate(aunitv, rotmat);
484                                        f.setUnitv(aunitv);
485
486                                        boolean doNotAdd = false;
487                                        if (params.reduceInitialFragments()) {
488                                                doNotAdd = FragmentJoiner.reduceFragments(
489                                                                fragments, f, distanceMatrix);
490
491                                        }
492                                        if (doNotAdd)
493                                                continue;
494
495                                        fragments.add(f);
496                                }
497                        }
498                }
499
500                notifyFragmentListeners(fragments);
501
502                FragmentPair[] fp = fragments
503                                .toArray(new FragmentPair[fragments.size()]);
504                setFragmentPairs(fp);
505
506                logger.debug(" got # fragment pairs: {}", fp.length);
507
508                logger.debug("step 2 - join fragments");
509
510                // step 2 combine them to possible models
511                FragmentJoiner joiner = new FragmentJoiner();
512
513                JointFragments[] frags;
514
515                if (params.isJoinFast()) {
516                        // apply the quick alignment procedure.
517                        // less quality in alignments, better for DB searches...
518                        frags = joiner.approach_ap3(ca1, ca2, fp, params);
519
520                        joiner.extendFragments(ca1, ca2, frags, params);
521
522                } else if (params.isJoinPlo()) {
523                        // this approach by StrComPy (peter lackner):
524                        frags = joiner.frag_pairwise_compat(fp, params.getAngleDiff(),
525                                        params.getFragCompat(), params.getMaxrefine());
526
527                } else {
528
529                        // my first implementation
530                        frags = joiner.approach_ap3(ca1, ca2, fp, params);
531                }
532
533                notifyJointFragments(frags);
534
535                logger.debug(" number joint fragments: ", frags.length);
536
537                logger.debug("step 3 - refine alignments");
538
539                List<AlternativeAlignment> aas = new ArrayList<AlternativeAlignment>();
540                for (int i = 0; i < frags.length; i++) {
541                        JointFragments f = frags[i];
542                        AlternativeAlignment a = new AlternativeAlignment();
543                        a.apairs_from_idxlst(f);
544                        a.setAltAligNumber(i + 1);
545                        a.setDistanceMatrix(distanceMatrix);
546
547                        try {
548                                if (params.getMaxIter() > 0) {
549
550                                        a.refine(params, ca1, ca2);
551                                } else {
552
553                                        a.finish(params, ca1, ca2);
554
555                                }
556                        } catch (StructureException e) {
557                                logger.error("Refinement of fragment {} failed", i, e);
558                        }
559                        a.calcScores(ca1, ca2);
560                        aas.add(a);
561                }
562
563                // sort the alternative alignments
564                Comparator<AlternativeAlignment> comp = new AltAligComparator();
565                Collections.sort(aas, comp);
566                Collections.reverse(aas);
567
568                alts = aas.toArray(new AlternativeAlignment[aas.size()]);
569                // do final numbering of alternative solutions
570                int aanbr = 0;
571                for (AlternativeAlignment a : alts) {
572                        aanbr++;
573                        a.setAltAligNumber(aanbr);
574                }
575
576                logger.debug("total calculation time: {} ms.",
577                                (System.currentTimeMillis() - timeStart));
578        }
579
580        private void notifyStartingAlignment(String name1, Atom[] ca1,
581                        String name2, Atom[] ca2) {
582                for (AlignmentProgressListener li : listeners) {
583                        li.startingAlignment(name1, ca1, name2, ca2);
584                }
585        }
586
587        private void notifyFragmentListeners(List<FragmentPair> fragments) {
588
589                for (AlignmentProgressListener li : listeners) {
590                        li.calculatedFragmentPairs(fragments);
591                }
592
593        }
594
595        private void notifyJointFragments(JointFragments[] fragments) {
596                for (AlignmentProgressListener li : listeners) {
597                        li.jointFragments(fragments);
598                }
599        }
600
601}