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

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

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