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

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

Building distance matrix between sequences
File size: 7.9 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 double 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                        double firstRowLeftScore = matrix[i-1][0].getScore();
104                        //for sequences of length 1
105                        double 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                        else {
142                        //End after we calculated final score
143                                return;
144                        }
145                       
146                       
147                        for (int j = 1; j < length2; j++) {
148                                double diagScore = matrix[i - 1][j - 1].getScore() + similarity(i, j);
149                                double upScore = matrix[i][j - 1].getScore() + submat.getGapPenalty();
150                                double leftScore = matrix[i - 1][j].getScore() + submat.getGapPenalty();
151
152                                matrix[i][j].setScore(Math.max(diagScore,Math.max(upScore, Math.max(leftScore,matrix[i][0].getScore()))));
153
154                                // find the directions that give the maximum scores.
155                                // Multiple directions are ignored TODO
156                                //True if we had a match
157                                if (diagScore == matrix[i][j].getScore()) {
158                                        matrix[i][j].setPrevious(matrix[i-1][j-1]);
159                                        matrix[i][j].setXvalue(input1[i-1]);
160                                        matrix[i][j].setYvalue(input2[j-1]);
161                                }
162                                //true if we took an event from sequence x and not from y
163                                if (leftScore == matrix[i][j].getScore()) {
164                                        matrix[i][j].setXvalue(input1[i-1]);
165                                        matrix[i][j].setYvalue(-1);
166                                        matrix[i][j].setPrevious(matrix[i-1][j]);
167                                }
168                                //true if we took an event from sequence y and not from x
169                                if (upScore == matrix[i][j].getScore()) {
170                                        matrix[i][j].setXvalue(-1);
171                                        matrix[i][j].setYvalue(input2[j-1]);
172                                        matrix[i][j].setPrevious(matrix[i][j-1]);
173                                }
174                                //true if we ended a matching region
175                                if (matrix[i][0].getScore() == matrix[i][j].getScore()) {
176                                        matrix[i][j].setPrevious(matrix[i][0]);
177                                        matrix[i][j].setXvalue(input1[i-1]);
178                                        matrix[i][j].setYvalue(-2);
179                                }
180                        }
181                       
182                        //Set the complete score cell
183                       
184                }
185        }
186
187        /**
188         * Get the maximum value in the score matrix.
189         */
190        public double getMaxScore() {
191                double maxScore = 0;
192
193                // skip the first row and column
194                for (int i = 1; i <= length1; i++) {
195                        for (int j = 1; j < length2; j++) {
196                                if (matrix[i][j].getScore() > maxScore) {
197                                        maxScore = matrix[i][j].getScore();
198                                }
199                        }
200                }
201
202                return maxScore;
203        }
204
205        /**
206         * Get the alignment score between the two input strings.
207         */
208        public double getAlignmentScore() {
209                return matrix[length1+1][0].getScore();
210        }
211
212       
213       
214        private int[] traceback(int i, int j) {
215               
216                       
217                return null;
218        }
219       
220        public List<Match> traceback() {
221                MatrixEntry tmp = matrix[length1+1][0];
222                String aligned1 = "";
223                String aligned2 = "";
224                int count = 0;
225                do
226                {
227                        String append1="";
228                        String append2="";
229                                       
230                        if(tmp.getXvalue() == -1) {
231                                append1 = "  ___";
232                        }
233                        else if(tmp.getXvalue() == -2) {
234                                append1 = "  ...";
235                        }
236                        else {
237                                append1 = String.format("%5d", tmp.getXvalue());
238                        }
239
240                        if(tmp.getYvalue() == -1) {
241                                append2 = "  ___";
242                        }
243                        else if(tmp.getYvalue() == -2) {
244                                append2 = "  ...";
245                        }
246                        else {
247                                append2 = String.format("%5d", tmp.getYvalue());
248                        }
249                        if(count != 0)
250                        {
251                                aligned1 = append1 + aligned1;
252                                aligned2 = append2 + aligned2;
253                        }
254                       
255                        tmp = tmp.getPrevious();
256                        count++;
257                       
258                } while(tmp != null);
259                System.out.println(aligned1);
260                System.out.println(aligned2);
261                return null;
262        }
263       
264
265       
266        /**
267         * print the dynmaic programming matrix
268         */
269        public void printDPMatrix() {
270                System.out.print("          ");
271                for (int i = 1; i <= length1; i++)
272                        System.out.format("%5d", input1[i - 1]);
273                System.out.println();
274                for (int j = 0; j < length2; j++) {
275                        if (j > 0)
276                                System.out.format("%5d ",input2[j - 1]);
277                        else{
278                                System.out.print("      ");
279                        }
280                        for (int i = 0; i <= length1 + 1; i++) {
281                                if((i<length1+1) || (i==length1+1 && j==0))     {
282                                        System.out.format("%4.1f ",matrix[i][j].getScore());
283                                }
284                               
285                        }                       
286                        System.out.println();
287                }
288        }
289
290        /**
291         * Return a set of Matches identified in Dynamic programming matrix. A match
292         * is a pair of subsequences whose score is higher than the preset
293         * scoreThreshold
294         **/
295        public List<Match> getMatches() {
296                ArrayList<Match> matchList = new ArrayList();
297                int fA = 0, fB = 0;
298                // skip the first row and column, find the next maxScore after
299                // prevmaxScore
300                for (int i = 1; i <= length1; i++) {
301                        for (int j = 1; j <= length2; j++) {
302                                if (matrix[i][j].getScore() > scoreThreshold
303                                                && matrix[i][j].getScore() > matrix[i - 1][j - 1].getScore()
304                                                && matrix[i][j].getScore() > matrix[i - 1][j].getScore()
305                                                && matrix[i][j].getScore() > matrix[i][j - 1].getScore()) {
306                                        if (i == length1 || j == length2
307                                                        || matrix[i][j].getScore() > matrix[i + 1][j + 1].getScore()) {
308                                                // should be lesser than prev maxScore
309                                                fA = i;
310                                                fB = j;
311                                                int[] f = traceback(fA, fB); // sets the x, y to
312                                                                                                                // startAlignment
313                                                                                                                // coordinates
314                                                System.out.println(f[0] + " " + i + " " + f[1] + " "
315                                                                + j + " " + matrix[i][j].getScore());
316                                                // TODO Add matches to matchList
317                                        }
318                                }
319                        }
320                }
321                return matchList;
322        }
323
324}
Note: See TracBrowser for help on using the repository browser.