source: trunk/autoquest-plugin-alignment/src/main/java/de/ugoe/cs/autoquest/plugin/alignment/SmithWaterman.java @ 1314

Last change on this file since 1314 was 1314, checked in by rkrimmel, 11 years ago

Hopefully resolved svn issues
File size: 8.2 KB
Line 
1package de.ugoe.cs.autoquest.plugin.alignment;
2
3import java.util.ArrayList;
4import java.util.List;
5
6public class SmithWaterman implements Alignment {
7
8        /**
9         * The first input
10         */
11        private int[] input1;
12
13        /**
14         * The second input String
15         */
16        private int[] input2;
17
18        /**
19         * The lengths of the input
20         */
21        private int length1, length2;
22
23        /**
24         * The score matrix. The true scores should be divided by the normalization
25         * factor.
26         */
27        private double[][] score;
28
29        /**
30         * Score threshold
31         */
32        private int scoreThreshold;;
33
34        /**
35         * The similarity function constants.
36         */
37        // static final int MATCH_SCORE = 10;
38        // static final int MISMATCH_SCORE = -8;
39        // static final int INDEL_SCORE = -9;
40
41        /**
42         * Constants of directions. Multiple directions are stored by bits. The zero
43         * direction is the starting point.
44         */
45        static final int DR_LEFT = 1; // 0001
46        static final int DR_UP = 2; // 0010
47        static final int DR_DIAG = 4; // 0100
48        static final int DR_ZERO = 8; // 1000
49
50        /**
51         * The directions pointing to the cells that give the maximum score at the
52         * current cell. The first index is the column index. The second index is
53         * the row index.
54         */
55        private int[][] prevCells;
56
57        /**
58         * Substitution matrix to calculate scores
59         */
60        private SubstitutionMatrix submat;
61
62        public SmithWaterman(int[] input1, int[] input2, SubstitutionMatrix submat) {
63                this.input1 = input1;
64                this.input2 = input2;
65                length1 = input1.length;
66                length2 = input2.length;
67                this.submat = submat;
68
69                System.out.println("Starting SmithWaterman algorithm with a "
70                                + submat.getClass() + " Substitution Matrix: " + submat.info());
71                scoreThreshold = 20;
72                score = new double[length1 + 1][length2 + 1];
73                prevCells = new int[length1 + 1][length2 + 1];
74
75                buildMatrix();
76        }
77
78        /**
79         * Compute the similarity score of substitution The position of the first
80         * character is 1. A position of 0 represents a gap.
81         *
82         * @param i
83         *            Position of the character in str1
84         * @param j
85         *            Position of the character in str2
86         * @return Cost of substitution of the character in str1 by the one in str2
87         */
88        private double similarity(int i, int j) {
89                if (i == 0 || j == 0) {
90                        // it's a gap (indel)
91                        return submat.getGapPenalty();
92                }
93                // System.out.println("Diag letters: " + input1[i-1] + " " +
94                // input2[j-1]);
95                // return (input1[i - 1] == input2[j - 1]) ? MATCH_SCORE :
96                // MISMATCH_SCORE;
97                return submat.getDistance(i - 1, j - 1);
98        }
99
100        /**
101         * Build the score matrix using dynamic programming. Note: The indel scores
102         * must be negative. Otherwise, the part handling the first row and column
103         * has to be modified.
104         */
105        private void buildMatrix() {
106                if (submat.getGapPenalty() >= 0) {
107                        throw new Error("Indel score must be negative");
108                }
109
110                int i; // length of prefix substring of str1
111                int j; // length of prefix substring of str2
112
113                // base case
114                score[0][0] = 0;
115                prevCells[0][0] = DR_ZERO; // starting point
116
117                // the first row
118                for (i = 1; i <= length1; i++) {
119                        score[i][0] = 0;
120                        prevCells[i][0] = DR_ZERO;
121                }
122
123                // the first column
124                for (j = 1; j <= length2; j++) {
125                        score[0][j] = 0;
126                        prevCells[0][j] = DR_ZERO;
127                }
128
129                // the rest of the matrix
130                for (i = 1; i <= length1; i++) {
131                        for (j = 1; j <= length2; j++) {
132                                double diagScore = score[i - 1][j - 1] + similarity(i, j);
133                                double upScore = score[i][j - 1] + similarity(0, j);
134                                double leftScore = score[i - 1][j] + similarity(i, 0);
135
136                                score[i][j] = Math.max(diagScore,
137                                                Math.max(upScore, Math.max(leftScore, 0)));
138                                prevCells[i][j] = 0;
139
140                                // find the directions that give the maximum scores.
141                                // the bitwise OR operator is used to record multiple
142                                // directions.
143                                if (diagScore == score[i][j]) {
144                                        prevCells[i][j] |= DR_DIAG;
145                                }
146                                if (leftScore == score[i][j]) {
147                                        prevCells[i][j] |= DR_LEFT;
148                                }
149                                if (upScore == score[i][j]) {
150                                        prevCells[i][j] |= DR_UP;
151                                }
152                                if (0 == score[i][j]) {
153                                        prevCells[i][j] |= DR_ZERO;
154                                }
155                        }
156                }
157        }
158
159        /**
160         * Get the maximum value in the score matrix.
161         */
162        private double getMaxScore() {
163                double maxScore = 0;
164
165                // skip the first row and column
166                for (int i = 1; i <= length1; i++) {
167                        for (int j = 1; j <= length2; j++) {
168                                if (score[i][j] > maxScore) {
169                                        maxScore = score[i][j];
170                                }
171                        }
172                }
173
174                return maxScore;
175        }
176
177        /**
178         * Get the alignment score between the two input strings.
179         */
180        public double getAlignmentScore() {
181                return getMaxScore();
182        }
183
184        /**
185         * TODO: Iterative Version!!! Output the local alignments ending in the (i,
186         * j) cell. aligned1 and aligned2 are suffixes of final aligned strings
187         * found in backtracking before calling this function. Note: the strings are
188         * replicated at each recursive call. Use buffers or stacks to improve
189         * efficiency.
190         */
191        private void printAlignments(int i, int j, String aligned1, String aligned2) {
192                // we've reached the starting point, so print the alignments
193
194                if ((prevCells[i][j] & DR_ZERO) > 0) {
195                        System.out.println(aligned1);
196                        System.out.println(aligned2);
197                        System.out.println();
198
199                        // Note: we could check other directions for longer alignments
200                        // with the same score. we don't do it here.
201                        return;
202                }
203
204                // find out which directions to backtrack
205                if ((prevCells[i][j] & DR_LEFT) > 0) {
206                        printAlignments(i - 1, j, input1[i - 1] + aligned1, "_ " + aligned2);
207                }
208                if ((prevCells[i][j] & DR_UP) > 0) {
209                        printAlignments(i, j - 1, "_ " + aligned1, input2[j - 1] + aligned2);
210                }
211                if ((prevCells[i][j] & DR_DIAG) > 0) {
212                        printAlignments(i - 1, j - 1, input1[i - 1] + " " + aligned1,
213                                        input2[j - 1] + " " + aligned2);
214                }
215        }
216
217        /**
218         * given the bottom right corner point trace back the top left conrner. at
219         * entry: i, j hold bottom right (end of Aligment coords) at return: hold
220         * top left (start of Alignment coords)
221         */
222        private int[] traceback(int i, int j) {
223
224                // find out which directions to backtrack
225                while (true) {
226                        if ((prevCells[i][j] & DR_LEFT) > 0) {
227                                if (score[i - 1][j] > 0)
228                                        i--;
229                                else
230                                        break;
231                        }
232                        if ((prevCells[i][j] & DR_UP) > 0) {
233                                // return traceback(i, j-1);
234                                if (score[i][j - 1] > 0)
235                                        j--;
236                                else
237                                        break;
238                        }
239                        if ((prevCells[i][j] & DR_DIAG) > 0) {
240                                // return traceback(i-1, j-1);
241                                if (score[i - 1][j - 1] > 0) {
242                                        i--;
243                                        j--;
244                                } else
245                                        break;
246                        }
247                }
248                int[] m = { i, j };
249                return m;
250        }
251
252        /**
253         * Output the local alignments with the maximum score.
254         */
255        public void printAlignments() {
256                // find the cell with the maximum score
257                double maxScore = getMaxScore();
258
259                /*
260                 * for (int i = 0; i < matches.length; i++) {
261                 * System.out.println("Match #" + i + ":" + matches.get(i)); }
262                 */
263
264                // skip the first row and column
265                for (int i = 1; i <= length1; i++) {
266                        for (int j = 1; j <= length2; j++) {
267                                if (score[i][j] == maxScore) {
268                                        printAlignments(i, j, "", "");
269                                }
270                        }
271                }
272                // Note: empty alignments are not printed.
273        }
274
275        /**
276         * print the dynmaic programming matrix
277         */
278        public void printDPMatrix() {
279                System.out.print("   ");
280                for (int j = 1; j <= length2; j++)
281                        System.out.print("   " + input2[j - 1]);
282                System.out.println();
283                for (int i = 0; i <= length1; i++) {
284                        if (i > 0)
285                                System.out.print(input1[i - 1] + " ");
286                        else
287                                System.out.print("  ");
288                        for (int j = 0; j <= length2; j++) {
289                                System.out.print(score[i][j] + " ");
290                        }
291                        System.out.println();
292                }
293        }
294
295        /**
296         * Return a set of Matches identified in Dynamic programming matrix. A match
297         * is a pair of subsequences whose score is higher than the preset
298         * scoreThreshold
299         **/
300        public List<Match> getMatches() {
301                ArrayList<Match> matchList = new ArrayList();
302                int fA = 0, fB = 0;
303                // skip the first row and column, find the next maxScore after
304                // prevmaxScore
305                for (int i = 1; i <= length1; i++) {
306                        for (int j = 1; j <= length2; j++) {
307                                if (score[i][j] > scoreThreshold
308                                                && score[i][j] > score[i - 1][j - 1]
309                                                && score[i][j] > score[i - 1][j]
310                                                && score[i][j] > score[i][j - 1]) {
311                                        if (i == length1 || j == length2
312                                                        || score[i][j] > score[i + 1][j + 1]) {
313                                                // should be lesser than prev maxScore
314                                                fA = i;
315                                                fB = j;
316                                                int[] f = traceback(fA, fB); // sets the x, y to
317                                                                                                                // startAlignment
318                                                                                                                // coordinates
319                                                System.out.println(f[0] + " " + i + " " + f[1] + " "
320                                                                + j + " " + score[i][j]);
321                                                // TODO Add matches to matchList
322                                        }
323                                }
324                        }
325                }
326                return matchList;
327        }
328
329}
Note: See TracBrowser for help on using the repository browser.