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 August 13, 2010
021 * Author: Mark Chapman
022 */
023
024package org.biojava.nbio.alignment.routines;
025
026import org.biojava.nbio.core.alignment.template.AlignedSequence.Step;
027
028import java.util.ArrayList;
029import java.util.Collections;
030import java.util.Comparator;
031import java.util.List;
032//import org.slf4j.Logger;
033//import org.slf4j.LoggerFactory;
034
035
036/**
037 * Static utility to construct alignment routines from a common library of methods.
038 *
039 * @author Mark Chapman
040 * @author Daniel Cameron
041 */
042public class AlignerHelper {
043
044        //private static final Logger logger = LoggerFactory.getLogger(AlignerHelper.class);
045
046        // types
047
048        /**
049         * Define a traceback pointer for the three edit operations: substitution (match/replacement of a query compound
050         * with a target compound), deletion (removal of a query compound leaving a gap in the target sequence), and
051         * insertion (addition of a target compound opening a gap in the query sequence).
052         */
053        public enum Last {
054                SUBSTITUTION,
055                DELETION,
056                INSERTION
057        }
058
059        /**
060         * Defines a 'cut' row for divide-and-conquer alignment in which a new anchor is found.
061         */
062        public static class Cut {
063
064                private int queryIndex;
065                private int[][] targetIndices, tiLast, ti1, ti2;
066
067                public Cut(int queryIndex, int[] dim) {
068                        this.queryIndex = queryIndex;
069                        targetIndices = ti1 = new int[dim[1]][dim[2]];
070                        ti2 = new int[dim[1]][dim[2]];
071                }
072
073                public int getQueryIndex() {
074                        return queryIndex;
075                }
076
077                public int getTargetIndex(int z) {
078                        return targetIndices[targetIndices.length - 1][z];
079                }
080
081                public void update(int x, Subproblem subproblem, Last[][] pointers) {
082                        if (pointers[subproblem.getTargetStartIndex()].length == 1) {
083                                if (queryIndex == x - 1) {
084                                        updateLinearInitial(subproblem, pointers);
085                                } else if (queryIndex < x) {
086                                        updateLinearAdvance(subproblem, pointers);
087                                }
088                        } else {
089                                if (queryIndex == x - 1) {
090                                        updateInitial(subproblem, pointers);
091                                } else if (queryIndex < x) {
092                                        updateAdvance(subproblem, pointers);
093                                }
094                        }
095                }
096
097                private void updateAdvance(Subproblem subproblem, Last[][] pointers) {
098                        tiLast = targetIndices;
099                        targetIndices = (targetIndices == ti2) ? ti1 : ti2;
100                        for (int y = subproblem.getTargetStartIndex(); y <= subproblem.getTargetEndIndex(); y++) {
101                                if (pointers[y][0] != null) {
102                                        targetIndices[y][0] = tiLast[y - 1][pointers[y][0].ordinal()];
103                                }
104                                if (pointers[y][1] != null) {
105                                        targetIndices[y][1] = tiLast[y][pointers[y][1].ordinal()];
106                                }
107                                if (pointers[y][2] != null) {
108                                        targetIndices[y][2] = targetIndices[y - 1][pointers[y][2].ordinal()];
109                                }
110                        }
111                }
112
113                private void updateInitial(Subproblem subproblem, Last[][] pointers) {
114                        for (int y = subproblem.getTargetStartIndex(); y <= subproblem.getTargetEndIndex(); y++) {
115                                if (pointers[y][0] != null) {
116                                        targetIndices[y][0] = y - 1;
117                                }
118                                if (pointers[y][1] != null) {
119                                        targetIndices[y][1] = y;
120                                }
121                                if (pointers[y][2] != null) {
122                                        targetIndices[y][2] = targetIndices[y - 1][2];
123                                }
124                        }
125                }
126
127                private void updateLinearAdvance(Subproblem subproblem, Last[][] pointers) {
128                        tiLast = targetIndices;
129                        targetIndices = (targetIndices == ti2) ? ti1 : ti2;
130                        for (int y = subproblem.getTargetStartIndex(); y <= subproblem.getTargetEndIndex(); y++) {
131                                switch (pointers[y][0]) {
132                                case DELETION:
133                                        targetIndices[y][0] = tiLast[y][0];
134                                        break;
135                                case SUBSTITUTION:
136                                        targetIndices[y][0] = tiLast[y - 1][0];
137                                        break;
138                                case INSERTION:
139                                        targetIndices[y][0] = targetIndices[y - 1][0];
140                                }
141                        }
142                }
143
144                private void updateLinearInitial(Subproblem subproblem, Last[][] pointers) {
145                        for (int y = subproblem.getTargetStartIndex(); y <= subproblem.getTargetEndIndex(); y++) {
146                                if (pointers[y][0] != null) {
147                                        switch (pointers[y][0]) {
148                                        case DELETION:
149                                                targetIndices[y][0] = y;
150                                                break;
151                                        case SUBSTITUTION:
152                                                targetIndices[y][0] = y - 1;
153                                                break;
154                                        case INSERTION:
155                                                targetIndices[y][0] = targetIndices[y - 1][0];
156                                        }
157                                }
158                        }
159                }
160
161        }
162        public static int addAnchors(Cut[] cuts, int[] scores, boolean addScore, int[] anchors) {
163                int zMax = 0, subscore = scores[0];
164                for (int z = 1; z < scores.length; z++) {
165                        if (scores[z] > subscore) {
166                                zMax = z;
167                                subscore = scores[z];
168                        }
169                }
170                for (Cut c : cuts) {
171                        anchors[c.getQueryIndex()] = c.getTargetIndex(zMax);
172                }
173                return addScore ? (int) subscore : 0;
174        }
175
176        public static Cut[] getCuts(int k, Subproblem subproblem, int[] dim, boolean anchor0) {
177                Cut[] cuts;
178                int m = subproblem.getQueryEndIndex() - subproblem.getQueryStartIndex() - (anchor0 ? 1 : 0);
179                if (k < m) {
180                        cuts = new Cut[k];
181                        int firstCutIndex = subproblem.getQueryStartIndex() + (anchor0 ? 1 : 0);
182                        for (int i = 0; i < k; i++) {
183                                cuts[i] = new Cut(firstCutIndex + i * (m - 1) / (k - 1), dim);
184                        }
185                } else {
186                        cuts = new Cut[m];
187                        for (int i = 0, x = subproblem.getQueryStartIndex() + (anchor0 ? 1 : 0); i < m; i++, x++) {
188                                cuts[i] = new Cut(x, dim);
189                        }
190                }
191                return cuts;
192        }
193        /**
194         * Compounds in query and target sequences that must align
195         * @author Daniel Cameron
196         */
197        public static class Anchor {
198                public int getQueryIndex() {
199                        return queryIndex;
200                }
201                public int getTargetIndex() {
202                        return targetIndex;
203                }
204                private final int queryIndex;
205                private final int targetIndex;
206                public Anchor(int queryIndex, int targetIndex) {
207                        this.queryIndex = queryIndex;
208                        this.targetIndex = targetIndex;
209                }
210                public static class QueryIndexComparator implements Comparator<Anchor> {
211                        @Override
212                        public int compare(Anchor o1, Anchor o2) {
213                                return o1.getQueryIndex() - o2.getQueryIndex();
214                        }
215                }
216        }
217        /**
218         * Alignment subproblem. The bounds of the subproblem are the
219         * indicies representing the inclusive bounds of the dynamic programming
220         * alignment problem.
221         * @author Daniel Cameron
222         */
223        public static class Subproblem {
224                public int getTargetStartIndex() {
225                        return targetStartIndex;
226                }
227                public int getQueryEndIndex() {
228                        return queryEndIndex;
229                }
230                public int getTargetEndIndex() {
231                        return targetEndIndex;
232                }
233                public int getQueryStartIndex() {
234                        return queryStartIndex;
235                }
236                /**
237                 * Indicates whether the start query and start target index compounds
238                 * are anchored to each other
239                 * @return true if the compounds are anchored in the alignment, false otherwise
240                 */
241                public boolean isStartAnchored() {
242                        return isAnchored;
243                }
244                private int queryStartIndex; // [0]
245                private int targetStartIndex; // [1]
246                private int queryEndIndex; // [2]
247                private int targetEndIndex; // [3]
248                private boolean isAnchored;
249                public Subproblem(int queryStartIndex, int targetStartIndex, int queryEndIndex, int targetEndIndex) {
250                        this(queryStartIndex, targetStartIndex, queryEndIndex, targetEndIndex, false);
251                }
252                public Subproblem(int queryStartIndex, int targetStartIndex, int queryEndIndex, int targetEndIndex, boolean isAnchored) {
253                        this.queryStartIndex = queryStartIndex;
254                        this.targetStartIndex = targetStartIndex;
255                        this.queryEndIndex = queryEndIndex;
256                        this.targetEndIndex = targetEndIndex;
257                        this.isAnchored = isAnchored;
258                }
259                /**
260                 * Convert a list of anchors into a subproblem list.
261                 * @param anchors anchored read pairs
262                 * @param querySequenceLength length of query sequence
263                 * @param targetSequenceLength length of target sequence
264                 * @return list alignment subproblems
265                 */
266                public static List<Subproblem> getSubproblems(List<Anchor> anchors, int querySequenceLength, int targetSequenceLength) {
267                        Collections.sort(anchors, new Anchor.QueryIndexComparator());
268                        List<Subproblem> list = new ArrayList<Subproblem>();
269                        Anchor last = new Anchor(-1, -1); // sentinal anchor
270                        boolean isAnchored = false;
271                        for (int i = 0; i < anchors.size(); i++) {
272                                if (anchors.get(i).targetIndex <= last.targetIndex ||
273                                        anchors.get(i).queryIndex <= last.queryIndex) {
274                                        throw new IllegalArgumentException("Anchor set must allow at least one possible alignment.");
275                                }
276                                list.add(new Subproblem(
277                                                last.queryIndex + 1,
278                                                last.targetIndex + 1,
279                                                anchors.get(i).queryIndex,
280                                                anchors.get(i).targetIndex,
281                                                isAnchored));
282                                last = anchors.get(i);
283                                isAnchored = true;
284                        }
285                        list.add(new Subproblem(
286                                        last.queryIndex + 1,
287                                        last.targetIndex + 1,
288                                        querySequenceLength,
289                                        targetSequenceLength,
290                                        isAnchored));
291                        return list;
292                }
293        }
294        // updates cut rows given the latest row of traceback pointers
295        public static void setCuts(int x, Subproblem subproblem, Last[][] pointers, Cut[]cuts) {
296                for (Cut c : cuts) {
297                        c.update(x, subproblem, pointers);
298                }
299        }
300        /**
301         * Calculate the optimal alignment score for the given sequence positions with an affine or constant gap penalty
302         * @param x position in query
303         * @param y position in target
304         * @param gop gap opening penalty
305         * @param gep gap extension penalty
306         * @param sub compound match score
307         * @param scores dynamic programming score matrix to fill at the given position
308         * @return traceback direction for substitution, deletion and insertion
309         */
310        public static Last[] setScorePoint(int x, int y, int gop, int gep, int sub, int[][][] scores) {
311                Last[] pointers = new Last[3];
312
313                // substitution
314                if (scores[x - 1][y - 1][1] >= scores[x - 1][y - 1][0] && scores[x - 1][y - 1][1] >= scores[x - 1][y - 1][2]) {
315                        scores[x][y][0] = scores[x - 1][y - 1][1] + sub;
316                        pointers[0] = Last.DELETION;
317                } else if (scores[x - 1][y - 1][0] >= scores[x - 1][y - 1][2]) {
318                        scores[x][y][0] = scores[x - 1][y - 1][0] + sub;
319                        pointers[0] = Last.SUBSTITUTION;
320                } else {
321                        scores[x][y][0] = scores[x - 1][y - 1][2] + sub;
322                        pointers[0] = Last.INSERTION;
323                }
324
325                // deletion
326                if (scores[x - 1][y][1] >= scores[x - 1][y][0] + gop) {
327                        scores[x][y][1] = scores[x - 1][y][1] + gep;
328                        pointers[1] = Last.DELETION;
329                } else {
330                        scores[x][y][1] = scores[x - 1][y][0] + gop + gep;
331                        pointers[1] = Last.SUBSTITUTION;
332                }
333
334                // insertion
335                if (scores[x][y - 1][0] + gop >= scores[x][y - 1][2]) {
336                        scores[x][y][2] = scores[x][y - 1][0] + gop + gep;
337                        pointers[2] = Last.SUBSTITUTION;
338                } else {
339                        scores[x][y][2] = scores[x][y - 1][2] + gep;
340                        pointers[2] = Last.INSERTION;
341                }
342
343                return pointers;
344        }
345        /**
346         * Calculates the optimal alignment score for the given sequence positions and a linear gap penalty
347         * @param x position in query
348         * @param y position in target
349         * @param gep gap extension penalty
350         * @param sub compound match score
351         * @param scores dynamic programming score matrix to fill at the given position
352         * @return traceback directions for substitution, deletion and insertion respectively
353         */
354        public static Last setScorePoint(int x, int y, int gep, int sub, int[][][] scores) {
355                int d = scores[x - 1][y][0] + gep;
356                int i = scores[x][y - 1][0] + gep;
357                int s = scores[x - 1][y - 1][0] + sub;
358                if (d >= s && d >= i) {
359                        scores[x][y][0] = d;
360                        return Last.DELETION;
361                } else if (s >= i) {
362                        scores[x][y][0] = s;
363                        return Last.SUBSTITUTION;
364                } else {
365                        scores[x][y][0] = i;
366                        return Last.INSERTION;
367                }
368        }
369
370        /**
371         * Score global alignment for a given position in the query sequence
372         * @param x
373         * @param subproblem
374         * @param gop
375         * @param gep
376         * @param subs
377         * @param storing
378         * @param scores
379         * @return
380         */
381        public static Last[][] setScoreVector(int x, Subproblem subproblem, int gop, int gep, int[] subs, boolean storing,
382                        int[][][] scores) {
383                return setScoreVector(x, subproblem.getQueryStartIndex(), subproblem.getTargetStartIndex(), subproblem.getTargetEndIndex(), gop, gep, subs, storing, scores, subproblem.isStartAnchored());
384        }
385
386        /**
387         * Score global alignment for a given position in the query sequence
388         * @param x
389         * @param xb
390         * @param yb
391         * @param ye
392         * @param gop
393         * @param gep
394         * @param subs
395         * @param storing
396         * @param scores
397         * @param startAnchored
398         * @return
399         */
400        public static Last[][] setScoreVector(int x, int xb, int yb, int ye, int gop, int gep, int[] subs,
401                        boolean storing, int[][][] scores, boolean startAnchored) {
402                Last[][] pointers = new Last[ye + 1][];
403                int min = Integer.MIN_VALUE - gop - gep;
404                ensureScoringMatrixColumn(x, storing, scores);
405                if (x == xb) {
406                        scores[xb][yb][1] = scores[xb][yb][2] = gop;
407                        pointers[yb] = new Last[] {null, null, null};
408                        if (startAnchored) {
409                                assert (xb > 0 && yb > 0);
410                                int subproblemStartingScore = scores[xb - 1][yb - 1][0] + subs[yb];
411                                scores[xb][yb][0] = subproblemStartingScore;
412                                scores[xb][yb][1] = subproblemStartingScore + gop;
413                                scores[xb][yb][2] = subproblemStartingScore + gop;
414                                pointers[yb] = new Last[] {Last.SUBSTITUTION, Last.SUBSTITUTION, Last.SUBSTITUTION};
415                        }
416                        Last[] insertion = new Last[] { null, null, Last.INSERTION };
417                        for (int y = yb + 1; y <= ye; y++) {
418                                scores[xb][y][0] = scores[xb][y][1] = min;
419                                scores[xb][y][2] = scores[xb][y - 1][2] + gep;
420                                pointers[y] = insertion;
421                        }
422                } else {
423                        scores[x][yb][0] = scores[x][yb][2] = min;
424                        scores[x][yb][1] = scores[x - 1][yb][1] + gep;
425                        pointers[yb] = new Last[] { null, Last.DELETION, null };
426                        for (int y = yb + 1; y <= ye; y++) {
427                                pointers[y] = setScorePoint(x, y, gop, gep, subs[y], scores);
428                        }
429                }
430                return pointers;
431        }
432
433        /**
434         * Score global alignment for a given position in the query sequence for a linear gap penalty
435         * @param x
436         * @param subproblem
437         * @param gep
438         * @param subs
439         * @param storing
440         * @param scores
441         * @return
442         */
443        public static Last[][] setScoreVector(int x, Subproblem subproblem, int gep, int[] subs, boolean storing,
444                        int[][][] scores) {
445                return setScoreVector(x, subproblem.getQueryStartIndex(), subproblem.getTargetStartIndex(), subproblem.getTargetEndIndex(), gep, subs, storing, scores, subproblem.isStartAnchored());
446        }
447
448        /**
449         * Score global alignment for a given position in the query sequence for a linear gap penalty
450         * @param x
451         * @param xb
452         * @param yb
453         * @param ye
454         * @param gep
455         * @param subs
456         * @param storing
457         * @param scores
458         * @param startAnchored
459         * @return
460         */
461        public static Last[][] setScoreVector(int x, int xb, int yb, int ye, int gep, int[] subs, boolean storing,
462                        int[][][] scores, boolean startAnchored) {
463                Last[][] pointers = new Last[ye + 1][1];
464                ensureScoringMatrixColumn(x, storing, scores);
465                if (x == xb) {
466                        if (startAnchored) {
467                                assert (xb > 0 && yb > 0);
468                                scores[xb][yb][0] = scores[xb - 1][yb - 1][0] + subs[yb];
469                                pointers[yb][0] = Last.SUBSTITUTION;
470                        }
471                        for (int y = yb + 1; y <= ye; y++) {
472                                scores[xb][y][0] = scores[xb][y - 1][0] + gep;
473                                pointers[y][0] = Last.INSERTION;
474                        }
475                } else {
476                        scores[x][yb][0] = scores[x - 1][yb][0] + gep;
477                        pointers[yb][0] = Last.DELETION;
478                        for (int y = yb + 1; y <= ye; y++) {
479                                pointers[y][0] = setScorePoint(x, y, gep, subs[y], scores);
480                        }
481                }
482                return pointers;
483        }
484
485        /**
486         * Score local alignment for a given position in the query sequence
487         * @param x
488         * @param gop
489         * @param gep
490         * @param subs
491         * @param storing
492         * @param scores
493         * @param xyMax
494         * @param score
495         * @return
496         */
497        public static Last[][] setScoreVector(int x, int gop, int gep, int[] subs, boolean storing,
498                        int[][][] scores, int[] xyMax, int score) {
499                return setScoreVector(x, 0, 0, scores[0].length - 1, gop, gep, subs, storing, scores, xyMax, score);
500        }
501
502        /**
503         * Score local alignment for a given position in the query sequence
504         * @param x
505         * @param xb
506         * @param yb
507         * @param ye
508         * @param gop
509         * @param gep
510         * @param subs
511         * @param storing
512         * @param scores
513         * @param xyMax
514         * @param score
515         * @return
516         */
517        public static Last[][] setScoreVector(int x, int xb, int yb, int ye, int gop, int gep, int[] subs,
518                        boolean storing, int[][][] scores, int[] xyMax, int score) {
519                Last[][] pointers;
520                ensureScoringMatrixColumn(x, storing, scores);
521                if (x == xb) {
522                        pointers = new Last[ye + 1][scores[0][0].length];
523                } else {
524                        pointers = new Last[ye + 1][];
525                        pointers[0] = new Last[scores[0][0].length];
526                        for (int y = 1; y < scores[0].length; y++) {
527                                pointers[y] = setScorePoint(x, y, gop, gep, subs[y], scores);
528                                for (int z = 0; z < scores[0][0].length; z++) {
529                                        if (scores[x][y][z] <= 0) {
530                                                scores[x][y][z] = 0;
531                                                pointers[y][z] = null;
532                                        }
533                                }
534                                if (scores[x][y][0] > score) {
535                                        xyMax[0] = x;
536                                        xyMax[1] = y;
537                                        score = scores[x][y][0];
538                                }
539                        }
540                }
541                return pointers;
542        }
543
544        /**
545         * Score local alignment for a given position in the query sequence for a linear gap penalty
546         * @param x
547         * @param gep
548         * @param subs
549         * @param storing
550         * @param scores
551         * @param xyMax
552         * @param score
553         * @return
554         */
555        public static Last[][] setScoreVector(int x, int gep, int[] subs, boolean storing, int[][][] scores,
556                        int[] xyMax, int score) {
557                return setScoreVector(x, 0, 0, scores[0].length - 1, gep, subs, storing, scores, xyMax, score);
558        }
559
560        /**
561         * Score local alignment for a given position in the query sequence for a linear gap penalty
562         * @param x
563         * @param xb
564         * @param yb
565         * @param ye
566         * @param gep
567         * @param subs
568         * @param storing
569         * @param scores
570         * @param xyMax
571         * @param score
572         * @return
573         */
574        public static Last[][] setScoreVector(int x, int xb, int yb, int ye, int gep, int[] subs, boolean storing,
575                        int[][][] scores, int[] xyMax, int score) {
576                Last[][] pointers;
577                ensureScoringMatrixColumn(x, storing, scores);
578                if (x == xb) {
579                        pointers = new Last[ye + 1][1];
580                } else {
581                        pointers = new Last[ye + 1][];
582                        pointers[0] = new Last[1];
583                        for (int y = 1; y < scores[x].length; y++) {
584                                pointers[y][0] = setScorePoint(x, y, gep, subs[y], scores);
585                                if (scores[x][y][0] <= 0) {
586                                        scores[x][y][0] = 0;
587                                        pointers[y][0] = null;
588                                } else if (scores[x][y][0] > score) {
589                                        xyMax[0] = x;
590                                        xyMax[1] = y;
591                                        score = scores[x][y][0];
592                                }
593                        }
594                }
595                return pointers;
596        }
597
598        private static void ensureScoringMatrixColumn(int x, boolean storingFullMatrix, int[][][] scores) {
599                if (!storingFullMatrix && x > 1) {
600                        scores[x] = scores[x - 2];
601                }
602        }
603
604        /**
605         * Find alignment path through traceback matrix
606         * @param traceback
607         * @param local
608         * @param xyMax
609         * @param last
610         * @param sx
611         * @param sy
612         * @return
613         */
614        public static int[] setSteps(Last[][][] traceback, boolean local, int[] xyMax, Last last, List<Step> sx,
615                        List<Step> sy) {
616                int x = xyMax[0], y = xyMax[1];
617                boolean linear = (traceback[x][y].length == 1);
618                while (local ? (linear ? last : traceback[x][y][last.ordinal()]) != null : x > 0 || y > 0) {
619                        switch (last) {
620                        case DELETION:
621                                sx.add(Step.COMPOUND);
622                                sy.add(Step.GAP);
623                                last = linear ? traceback[--x][y][0] : traceback[x--][y][1];
624                                break;
625                        case SUBSTITUTION:
626                                sx.add(Step.COMPOUND);
627                                sy.add(Step.COMPOUND);
628                                last = linear ? traceback[--x][--y][0] : traceback[x--][y--][0];
629                                break;
630                        case INSERTION:
631                                sx.add(Step.GAP);
632                                sy.add(Step.COMPOUND);
633                                last = linear ? traceback[x][--y][0] : traceback[x][y--][2];
634                        }
635                }
636                Collections.reverse(sx);
637                Collections.reverse(sy);
638                return new int[] {x, y};
639        }
640
641        /**
642         * Find global alignment path through traceback matrix
643         * @param traceback
644         * @param scores
645         * @param sx
646         * @param sy
647         * @return
648         */
649        public static int[] setSteps(Last[][][] traceback, int[][][] scores, List<Step> sx, List<Step> sy) {
650                int xMax = scores.length - 1, yMax = scores[xMax].length - 1;
651                boolean linear = (traceback[xMax][yMax].length == 1);
652
653                Last last =
654
655                        linear ?
656                                traceback[xMax][yMax][0] :
657
658                                (scores[xMax][yMax][1] > scores[xMax][yMax][0] &&
659                                 scores[xMax][yMax][1] > scores[xMax][yMax][2] ) ?
660
661                                                Last.DELETION :
662                                                        (scores[xMax][yMax][0] > scores[xMax][yMax][2]) ?
663                                                                        Last.SUBSTITUTION :
664                                                                        Last.INSERTION;
665
666
667                return setSteps(traceback, false, new int[] {xMax, yMax}, last, sx, sy);
668        }
669
670        /**
671         * Find local alignment path through traceback matrix
672         * @param traceback
673         * @param xyMax
674         * @param sx
675         * @param sy
676         * @return
677         */
678        public static int[] setSteps(Last[][][] traceback, int[] xyMax, List<Step> sx, List<Step> sy) {
679                return setSteps(traceback, true, xyMax, Last.SUBSTITUTION, sx, sy);
680        }
681
682        public static String tracebackToString(Last[][][] traceback) {
683                StringBuilder sb = new StringBuilder();
684                for (int z = 0; z < 3; z++) {
685                        for (int i = 0; i < traceback.length; i++) {
686                                if (traceback[i] != null) {
687                                        for (int j = 0; j < traceback[i].length; j++) {
688                                                if (traceback[i][j] == null || z >= traceback[i][j].length || traceback[i][j][z] == null) {
689                                                        sb.append('.');
690                                                } else {
691                                                        switch (traceback[i][j][z]) {
692                                                        case DELETION:
693                                                                sb.append('^');
694                                                                break;
695                                                        case SUBSTITUTION:
696                                                                sb.append('\\');
697                                                                break;
698                                                        case INSERTION:
699                                                                sb.append('<');
700                                                                break;
701                                                        }
702                                                }
703                                        }
704                                }
705                                sb.append('\n');
706                        }
707                        sb.append("\n\n");
708                }
709                return sb.toString();
710        }
711}