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

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

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