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