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

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

Multiple Alignment first version
File size: 8.4 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 {
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        /**
211         * Get the alignment score between the two input strings.
212         */
213        public double getAlignmentScore() {
214                return matrix[length1+1][0].getScore();
215        }
216
217       
218       
219        public void traceback() {
220                MatrixEntry tmp = matrix[length1+1][0];
221                int aligned1[] = new int[length1+length2+2];
222                int aligned2[] = new int[length1+length2+2];
223                int count = 0;
224                do
225                {       
226                        if(count != 0)
227                        {
228                                if (length1+length2+2 == count) {       
229                                        Console.traceln(Level.WARNING, "Traceback longer than both sequences summed up!");
230                                        break;
231                                }
232                                aligned1[count] = tmp.getXvalue();
233                                aligned2[count] = tmp.getYvalue();
234                        }
235                       
236                        tmp = tmp.getPrevious();
237                        count++;
238               
239                } while(tmp != null);
240                count--;
241                //reverse order of the alignment
242                int reversed1[] = new int[count];
243                int reversed2[] = new int[count];
244               
245               
246                for(int i = count-1; i > 0; i--) {
247                        reversed1[reversed1.length-i]= aligned1[i];
248                        reversed2[reversed2.length-i]= aligned2[i];
249                }
250               
251                NumberSequence ns1 = new NumberSequence(reversed1.length);
252                NumberSequence ns2 = new NumberSequence(reversed2.length);
253                ns1.setSequence(reversed1);
254                ns2.setSequence(reversed2);
255               
256                alignment.add(ns1);
257                alignment.add(ns2);
258        }
259       
260        public void printAlignment() {
261                MatrixEntry tmp = matrix[length1+1][0];
262                String aligned1 = "";
263                String aligned2 = "";
264                int count = 0;
265                do
266                {
267                        String append1="";
268                        String append2="";
269                                       
270                        if(tmp.getXvalue() == Constants.GAP_SYMBOL) {
271                                append1 = "  ___";
272                        }
273                        else if(tmp.getXvalue() == Constants.UNMATCHED_SYMBOL) {
274                                append1 = "  ...";
275                        }
276                        else {
277                                append1 = String.format("%5d", tmp.getXvalue());
278                        }
279
280                        if(tmp.getYvalue() == Constants.GAP_SYMBOL) {
281                                append2 = "  ___";
282                        }
283                        else if(tmp.getYvalue() == Constants.UNMATCHED_SYMBOL) {
284                                append2 = "  ...";
285                        }
286                        else {
287                                append2 = String.format("%5d", tmp.getYvalue());
288                        }
289                        if(count != 0)
290                        {
291                                aligned1 = append1 + aligned1;
292                                aligned2 = append2 + aligned2;
293                        }
294                       
295                        tmp = tmp.getPrevious();
296                        count++;
297                       
298                } while(tmp != null);
299                System.out.println(aligned1);
300                System.out.println(aligned2);
301        }
302       
303
304       
305        /**
306         * print the dynmaic programming matrix
307         */
308        public void printDPMatrix() {
309                System.out.print("          ");
310                for (int i = 1; i <= length1; i++)
311                        System.out.format("%5d", input1[i - 1]);
312                System.out.println();
313                for (int j = 0; j < length2; j++) {
314                        if (j > 0)
315                                System.out.format("%5d ",input2[j - 1]);
316                        else{
317                                System.out.print("      ");
318                        }
319                        for (int i = 0; i <= length1 + 1; i++) {
320                                if((i<length1+1) || (i==length1+1 && j==0))     {
321                                        System.out.format("%4.1f ",matrix[i][j].getScore());
322                                }
323                               
324                        }                       
325                        System.out.println();
326                }
327        }
328
329
330        public ArrayList<NumberSequence> getAlignment() {
331                return alignment;
332        }
333
334        public void setAlignment(ArrayList<NumberSequence> alignment) {
335                this.alignment = alignment;
336        }
337
338}
Note: See TracBrowser for help on using the repository browser.