source: branches/ralph/src/main/java/de/ugoe/cs/autoquest/tasktrees/alignment/algorithms/SmithWaterman.java @ 1586

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

Adding simple smith waterman and changing alignment algorithm creation to factory pattern
File size: 7.2 KB
Line 
1package de.ugoe.cs.autoquest.tasktrees.alignment.algorithms;
2
3import java.util.ArrayList;
4import java.util.logging.Level;
5
6import de.ugoe.cs.autoquest.tasktrees.alignment.matrix.SubstitutionMatrix;
7import de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.Constants;
8import de.ugoe.cs.util.console.Console;
9
10public class SmithWaterman implements AlignmentAlgorithm {
11
12        /**
13         * The first input
14         */
15        private int[] input1;
16
17        /**
18         * The second input String
19         */
20        private int[] input2;
21
22        /**
23         * The lengths of the input
24         */
25        private int length1, length2;
26
27        /**
28         * The score matrix. The true scores should be divided by the normalization
29         * factor.
30         */
31        private MatrixEntry[][] matrix;
32
33        private ArrayList<NumberSequence> alignment;
34
35        private float scoreThreshold;
36
37        /**
38         * Substitution matrix to calculate scores
39         */
40        private SubstitutionMatrix submat;
41
42        public SmithWaterman(int[] input1, int[] input2, SubstitutionMatrix submat,
43                        float threshold) {
44                this.input1 = input1;
45                this.input2 = input2;
46                length1 = input1.length;
47                length2 = input2.length;
48                this.submat = submat;
49
50                // System.out.println("Starting SmithWaterman algorithm with a "
51                // + submat.getClass() + " Substitution Matrix: " +
52                // submat.getClass().getCanonicalName());
53                this.scoreThreshold = threshold;
54
55                matrix = new MatrixEntry[length1 + 1][length2 + 1];
56                alignment = new ArrayList<NumberSequence>();
57
58                for (int i = 0; i < length1+1; i++) {
59                        for (int j = 0; j < length2+1; j++) {
60                                matrix[i][j] = new MatrixEntry();
61                        }
62                }
63
64                buildMatrix();
65                traceback();
66        }
67
68        /**
69         * Compute the similarity score of substitution The position of the first
70         * character is 1. A position of 0 represents a gap.
71         *
72         * @param i
73         *            Position of the character in str1
74         * @param j
75         *            Position of the character in str2
76         * @return Cost of substitution of the character in str1 by the one in str2
77         */
78        private double similarity(int i, int j) {
79                return submat.getScore(input1[i - 1], input2[j - 1]);
80        }
81
82        /**
83         * Build the score matrix using dynamic programming.
84         */
85        private void buildMatrix() {
86                if (submat.getGapPenalty() >= 0) {
87                        throw new Error("Indel score must be negative");
88                }
89
90                // it's a gap
91                matrix[0][0].setScore(0);
92                matrix[0][0].setPrevious(null); // starting point
93
94                // the first column
95                for (int j = 1; j < length2+1; j++) {
96                        matrix[0][j].setScore(0);
97                        matrix[0][j].setPrevious(matrix[0][j - 1]);
98                        matrix[0][j].setYvalue(input2[j]);
99                        matrix[0][j].setXvalue(Constants.GAP_SYMBOL);
100                }
101                // the first row
102
103                for (int j = 1; j < length1+1; j++) {
104                        matrix[j][0].setScore(0);
105                        matrix[j][0].setPrevious(matrix[j - 1][0]);
106                        matrix[j][0].setXvalue(input1[j]);
107                        matrix[j][0].setYvalue(Constants.GAP_SYMBOL);
108                }
109
110                for (int i = 1; i < length1; i++) {
111
112                        for (int j = 1; j < length2; j++) {
113                                double diagScore = matrix[i - 1][j - 1].getScore()
114                                                + similarity(i, j);
115                                double upScore = matrix[i][j - 1].getScore()
116                                                + submat.getGapPenalty();
117                                double leftScore = matrix[i - 1][j].getScore()
118                                                + submat.getGapPenalty();
119
120                                matrix[i][j].setScore(Math.max(diagScore,
121                                                Math.max(upScore, Math.max(leftScore, 0))));
122
123                                // find the directions that give the maximum scores.
124                                // TODO: Multiple directions are ignored, we choose the first
125                                // maximum score
126                                // True if we had a match
127                                if (diagScore == matrix[i][j].getScore()) {
128                                        matrix[i][j].setPrevious(matrix[i - 1][j - 1]);
129                                        matrix[i][j].setXvalue(input1[i - 1]);
130                                        matrix[i][j].setYvalue(input2[j - 1]);
131                                }
132                                // true if we took an event from sequence x and not from y
133                                if (leftScore == matrix[i][j].getScore()) {
134                                        matrix[i][j].setXvalue(input1[i - 1]);
135                                        matrix[i][j].setYvalue(Constants.GAP_SYMBOL);
136                                        matrix[i][j].setPrevious(matrix[i - 1][j]);
137                                }
138                                // true if we took an event from sequence y and not from x
139                                if (upScore == matrix[i][j].getScore()) {
140                                        matrix[i][j].setXvalue(Constants.GAP_SYMBOL);
141                                        matrix[i][j].setYvalue(input2[j - 1]);
142                                        matrix[i][j].setPrevious(matrix[i][j - 1]);
143                                }
144                                if (0 == matrix[i][j].getScore()) {
145                                        matrix[i][j].setPrevious(matrix[0][0]);
146                                        matrix[i][j].setXvalue(-2);
147                                        matrix[i][j].setYvalue(-2);
148                                }
149                        }
150
151                        // Set the complete score cell
152
153                }
154        }
155
156        /**
157         * Get the maximum value in the score matrix.
158         */
159        public double getMaxScore() {
160                double maxScore = 0;
161
162                // skip the first row and column
163                for (int i = 1; i <= length1; i++) {
164                        for (int j = 1; j < length2; j++) {
165                                if (matrix[i][j].getScore() > maxScore) {
166                                        maxScore = matrix[i][j].getScore();
167                                }
168                        }
169                }
170
171                return maxScore;
172        }
173
174        /*
175         * (non-Javadoc)
176         *
177         * @see
178         * de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.AlignmentAlgorithm
179         * #getAlignmentScore()
180         */
181        @Override
182        public double getAlignmentScore() {
183                return getMaxScore();
184        }
185
186        public void traceback() {
187                MatrixEntry tmp = matrix[length1][length2];
188                int aligned1[] = new int[length1 + length2 + 2];
189                int aligned2[] = new int[length1 + length2 + 2];
190                int count = 0;
191                do {
192                        if (length1 + length2 + 2 == count) {
193                                Console.traceln(Level.WARNING,
194                                                "Traceback longer than both sequences summed up!");
195                                break;
196                        }
197                        aligned1[count] = tmp.getXvalue();
198                        aligned2[count] = tmp.getYvalue();
199
200                        tmp = tmp.getPrevious();
201                        count++;
202
203                } while (tmp != null);
204                count--;
205                // reverse order of the alignment
206                int reversed1[] = new int[count];
207                int reversed2[] = new int[count];
208
209                for (int i = count - 1; i > 0; i--) {
210                        reversed1[reversed1.length - i] = aligned1[i];
211                        reversed2[reversed2.length - i] = aligned2[i];
212                }
213
214                NumberSequence ns1 = new NumberSequence(reversed1.length);
215                NumberSequence ns2 = new NumberSequence(reversed2.length);
216                ns1.setSequence(reversed1);
217                ns2.setSequence(reversed2);
218
219                alignment.add(ns1);
220                alignment.add(ns2);
221        }
222
223        public void printAlignment() {
224                MatrixEntry tmp = matrix[length1 + 1][0];
225                String aligned1 = "";
226                String aligned2 = "";
227                int count = 0;
228                do {
229                        String append1 = "";
230                        String append2 = "";
231
232                        if (tmp.getXvalue() == Constants.GAP_SYMBOL) {
233                                append1 = "  ___";
234                        } else {
235                                append1 = String.format("%5d", tmp.getXvalue());
236                        }
237
238                        if (tmp.getYvalue() == Constants.GAP_SYMBOL) {
239                                append2 = "  ___";
240                        } else {
241                                append2 = String.format("%5d", tmp.getYvalue());
242                        }
243                        if (count != 0) {
244                                aligned1 = append1 + aligned1;
245                                aligned2 = append2 + aligned2;
246                        }
247
248                        tmp = tmp.getPrevious();
249                        count++;
250
251                } while (tmp != null);
252                System.out.println(aligned1);
253                System.out.println(aligned2);
254        }
255
256        /**
257         * print the dynmaic programming matrix
258         */
259        public void printDPMatrix() {
260                System.out.print("          ");
261                for (int i = 1; i <= length1; i++)
262                        System.out.format("%5d", input1[i - 1]);
263                System.out.println();
264                for (int j = 0; j < length2; j++) {
265                        if (j > 0)
266                                System.out.format("%5d ", input2[j - 1]);
267                        else {
268                                System.out.print("      ");
269                        }
270                        for (int i = 0; i < length1; i++) {
271                                        System.out.format("%4.1f ", matrix[i][j].getScore());
272                        }
273                        System.out.println();
274                }
275        }
276
277        /*
278         * (non-Javadoc)
279         *
280         * @see
281         * de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.AlignmentAlgorithm
282         * #getAlignment()
283         */
284        @Override
285        public ArrayList<NumberSequence> getAlignment() {
286                return alignment;
287        }
288
289        public void setAlignment(ArrayList<NumberSequence> alignment) {
290                this.alignment = alignment;
291        }
292
293}
Note: See TracBrowser for help on using the repository browser.