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

Last change on this file since 1323 was 1323, checked in by rkrimmel, 10 years ago

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