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

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

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