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