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

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

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