source: branches/ralph/src/main/java/de/ugoe/cs/autoquest/tasktrees/alignment/algorithms/SmithWatermanRepeated.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: 8.6 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 SmithWatermanRepeated 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
34        private ArrayList<NumberSequence> alignment;
35       
36        private float scoreThreshold;
37
38        /**
39         * Substitution matrix to calculate scores
40         */
41        private SubstitutionMatrix submat;
42
43        public SmithWatermanRepeated(int[] input1, int[] input2, SubstitutionMatrix submat,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: " + submat.getClass().getCanonicalName());
52                this.scoreThreshold = threshold;
53               
54                matrix = new MatrixEntry[length1+2][length2+1];
55                alignment = new ArrayList<NumberSequence>();
56               
57                for (int i = 0; i <= length1+1; i++) {
58                        for(int j = 0; j< length2; j++) {
59                                matrix[i][j] = new MatrixEntry();
60                        }
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; j++) {
96                        matrix[0][j].setScore(0);
97                        //We don't need to go back to [0][0] if we reached matrix[0][x], so just end here
98                        //matrix[0][j].setPrevious(matrix[0][j-1]);
99                        matrix[0][j].setPrevious(null);
100                }
101               
102               
103               
104                for (int i = 1; i <= length1 + 1; i++) {
105               
106                        // Formula for first row:
107                        // F(i,0) = max { F(i-1,0), F(i-1,j)-T  j=1,...,m
108                       
109                        double firstRowLeftScore = matrix[i-1][0].getScore();
110                        //for sequences of length 1
111                        double tempMax;
112                        int maxRowIndex;
113                        if(length2 == 1) {
114                                tempMax = matrix[i-1][0].getScore();
115                                maxRowIndex = 0;
116                        } else {
117                                tempMax = matrix[i-1][1].getScore();
118                                maxRowIndex = 1;
119                                //position of the maximal score of the previous row
120                               
121                                for(int j = 2; j < length2;j++) {
122                                        if(matrix[i-1][j].getScore() > tempMax) {
123                                                tempMax = matrix[i-1][j].getScore();
124                                                maxRowIndex = j;
125                                        }
126                                }
127
128                        }
129                                       
130                        tempMax -= scoreThreshold;
131                        matrix[i][0].setScore(Math.max(firstRowLeftScore, tempMax));
132                        if(tempMax ==matrix[i][0].getScore()){
133                                matrix[i][0].setPrevious(matrix[i-1][maxRowIndex]);
134                        }
135                       
136                        if(firstRowLeftScore == matrix[i][0].getScore()) {
137                                matrix[i][0].setPrevious(matrix[i-1][0]);
138                        }
139                       
140                        //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
141                        if(i<length1+1)
142                        {
143                                matrix[i][0].setXvalue(input1[i-1]);
144                                matrix[i][0].setYvalue(Constants.UNMATCHED_SYMBOL);
145                        }
146                        else {
147                        //End after we calculated final score
148                                return;
149                        }
150                       
151                       
152                        for (int j = 1; j < length2; j++) {
153                                double diagScore = matrix[i - 1][j - 1].getScore() + similarity(i, j);
154                                double upScore = matrix[i][j - 1].getScore() + submat.getGapPenalty();
155                                double leftScore = matrix[i - 1][j].getScore() + submat.getGapPenalty();
156
157                                matrix[i][j].setScore(Math.max(diagScore,Math.max(upScore, Math.max(leftScore,matrix[i][0].getScore()))));
158
159                                // find the directions that give the maximum scores.
160                                // TODO: Multiple directions are ignored, we choose the first maximum score
161                                //True if we had a match
162                                if (diagScore == matrix[i][j].getScore()) {
163                                        matrix[i][j].setPrevious(matrix[i-1][j-1]);
164                                        matrix[i][j].setXvalue(input1[i-1]);
165                                        matrix[i][j].setYvalue(input2[j-1]);
166                                }
167                                //true if we took an event from sequence x and not from y
168                                if (leftScore == matrix[i][j].getScore()) {
169                                        matrix[i][j].setXvalue(input1[i-1]);
170                                        matrix[i][j].setYvalue(Constants.GAP_SYMBOL);
171                                        matrix[i][j].setPrevious(matrix[i-1][j]);
172                                }
173                                //true if we took an event from sequence y and not from x
174                                if (upScore == matrix[i][j].getScore()) {
175                                        matrix[i][j].setXvalue(Constants.GAP_SYMBOL);
176                                        matrix[i][j].setYvalue(input2[j-1]);
177                                        matrix[i][j].setPrevious(matrix[i][j-1]);
178                                }
179                                //true if we ended a matching region
180                                if (matrix[i][0].getScore() == matrix[i][j].getScore()) {
181                                        matrix[i][j].setPrevious(matrix[i][0]);
182                                        matrix[i][j].setXvalue(input1[i-1]);
183                                        matrix[i][j].setYvalue(Constants.UNMATCHED_SYMBOL);
184                                }
185                        }
186                       
187                        //Set the complete score cell
188                       
189                }
190        }
191
192        /**
193         * Get the maximum value in the score matrix.
194         */
195        public double getMaxScore() {
196                double maxScore = 0;
197
198                // skip the first row and column
199                for (int i = 1; i <= length1; i++) {
200                        for (int j = 1; j < length2; j++) {
201                                if (matrix[i][j].getScore() > maxScore) {
202                                        maxScore = matrix[i][j].getScore();
203                                }
204                        }
205                }
206
207                return maxScore;
208        }
209
210        /* (non-Javadoc)
211         * @see de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.AlignmentAlgorithm#getAlignmentScore()
212         */
213        @Override
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        /* (non-Javadoc)
332         * @see de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.AlignmentAlgorithm#getAlignment()
333         */
334        @Override
335        public ArrayList<NumberSequence> getAlignment() {
336                return alignment;
337        }
338
339        public void setAlignment(ArrayList<NumberSequence> alignment) {
340                this.alignment = alignment;
341        }
342
343}
Note: See TracBrowser for help on using the repository browser.