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

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

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