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

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

Added parts of the PAL library, implemented UPGMA Algoritm for Feng Doolittle guide tree
File size: 8.9 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;
7
8import de.ugoe.cs.autoquest.tasktrees.alignment.matrix.SubstitutionMatrix;
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 List<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        }
66
67        /**
68         * Compute the similarity score of substitution The position of the first
69         * character is 1. A position of 0 represents a gap.
70         *
71         * @param i
72         *            Position of the character in str1
73         * @param j
74         *            Position of the character in str2
75         * @return Cost of substitution of the character in str1 by the one in str2
76         */
77        private double similarity(int i, int j) {
78                return submat.getDistance(input1[i - 1], input2[j - 1]);
79        }
80
81        /**
82         * Build the score matrix using dynamic programming.
83         */
84        private void buildMatrix() {
85                if (submat.getGapPenalty() >= 0) {
86                        throw new Error("Indel score must be negative");
87                }
88               
89                // it's a gap
90                matrix[0][0].setScore(0);
91                matrix[0][0].setPrevious(null); // starting point
92
93                // the first column
94                for (int j = 1; j < length2; j++) {
95                        matrix[0][j].setScore(0);
96                        matrix[0][j].setPrevious(matrix[0][j-1]);
97                       
98                }
99               
100               
101               
102                for (int i = 1; i <= length1 + 1; i++) {
103               
104                        // Formula for first row:
105                        // F(i,0) = max { F(i-1,0), F(i-1,j)-T  j=1,...,m
106                       
107                        double firstRowLeftScore = matrix[i-1][0].getScore();
108                        //for sequences of length 1
109                        double tempMax;
110                        int maxRowIndex;
111                        if(length2 == 1) {
112                                tempMax = matrix[i-1][0].getScore();
113                                maxRowIndex = 0;
114                        } else {
115                                tempMax = matrix[i-1][1].getScore();
116                                maxRowIndex = 1;
117                                //position of the maximal score of the previous row
118                               
119                                for(int j = 2; j < length2;j++) {
120                                        if(matrix[i-1][j].getScore() > tempMax) {
121                                                tempMax = matrix[i-1][j].getScore();
122                                                maxRowIndex = j;
123                                        }
124                                }
125
126                        }
127                       
128                                       
129                        tempMax -= scoreThreshold;
130                        matrix[i][0].setScore(Math.max(firstRowLeftScore, tempMax));
131                        if(tempMax ==matrix[i][0].getScore()){
132                                matrix[i][0].setPrevious(matrix[i-1][maxRowIndex]);
133                        }
134                       
135                        if(firstRowLeftScore == matrix[i][0].getScore()) {
136                                matrix[i][0].setPrevious(matrix[i-1][0]);
137                        }
138                       
139                        //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
140                        if(i<length1+1)
141                        {
142                                matrix[i][0].setXvalue(input1[i-1]);
143                                matrix[i][0].setYvalue(-2);
144                        }
145                        else {
146                        //End after we calculated final score
147                                return;
148                        }
149                       
150                       
151                        for (int j = 1; j < length2; j++) {
152                                double diagScore = matrix[i - 1][j - 1].getScore() + similarity(i, j);
153                                double upScore = matrix[i][j - 1].getScore() + submat.getGapPenalty();
154                                double leftScore = matrix[i - 1][j].getScore() + submat.getGapPenalty();
155
156                                matrix[i][j].setScore(Math.max(diagScore,Math.max(upScore, Math.max(leftScore,matrix[i][0].getScore()))));
157
158                                // find the directions that give the maximum scores.
159                                // Multiple directions are ignored TODO
160                                //True if we had a match
161                                if (diagScore == matrix[i][j].getScore()) {
162                                        matrix[i][j].setPrevious(matrix[i-1][j-1]);
163                                        matrix[i][j].setXvalue(input1[i-1]);
164                                        matrix[i][j].setYvalue(input2[j-1]);
165                                }
166                                //true if we took an event from sequence x and not from y
167                                if (leftScore == matrix[i][j].getScore()) {
168                                        matrix[i][j].setXvalue(input1[i-1]);
169                                        matrix[i][j].setYvalue(-1);
170                                        matrix[i][j].setPrevious(matrix[i-1][j]);
171                                }
172                                //true if we took an event from sequence y and not from x
173                                if (upScore == matrix[i][j].getScore()) {
174                                        matrix[i][j].setXvalue(-1);
175                                        matrix[i][j].setYvalue(input2[j-1]);
176                                        matrix[i][j].setPrevious(matrix[i][j-1]);
177                                }
178                                //true if we ended a matching region
179                                if (matrix[i][0].getScore() == matrix[i][j].getScore()) {
180                                        matrix[i][j].setPrevious(matrix[i][0]);
181                                        matrix[i][j].setXvalue(input1[i-1]);
182                                        matrix[i][j].setYvalue(-2);
183                                }
184                        }
185                       
186                        //Set the complete score cell
187                       
188                }
189        }
190
191        /**
192         * Get the maximum value in the score matrix.
193         */
194        public double getMaxScore() {
195                double maxScore = 0;
196
197                // skip the first row and column
198                for (int i = 1; i <= length1; i++) {
199                        for (int j = 1; j < length2; j++) {
200                                if (matrix[i][j].getScore() > maxScore) {
201                                        maxScore = matrix[i][j].getScore();
202                                }
203                        }
204                }
205
206                return maxScore;
207        }
208
209        /**
210         * Get the alignment score between the two input strings.
211         */
212        public double getAlignmentScore() {
213                return matrix[length1+1][0].getScore();
214        }
215
216       
217       
218        private int[] traceback(int i, int j) {
219               
220                       
221                return null;
222        }
223       
224       
225        public void traceback() {
226                MatrixEntry tmp = matrix[length1+1][0];
227               
228                int aligned1[] = new int[length1+length2];
229                int aligned2[] = new int[length1+length2];
230                int count = 0;
231                do
232                {       
233                        if(count != 0)
234                        {
235                                aligned1[count] = tmp.getXvalue();
236                                aligned2[count] = tmp.getYvalue();
237                        }
238                       
239                        tmp = tmp.getPrevious();
240                        count++;
241                       
242                } while(tmp != null);
243
244                //reverse order
245                int reversed1[] = new int[count];
246                int reversed2[] = new int[count];
247               
248                for(int i = count; count > 0; count ++) {
249                       
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() == -1) {
272                                append1 = "  ___";
273                        }
274                        else if(tmp.getXvalue() == -2) {
275                                append1 = "  ...";
276                        }
277                        else {
278                                append1 = String.format("%5d", tmp.getXvalue());
279                        }
280
281                        if(tmp.getYvalue() == -1) {
282                                append2 = "  ___";
283                        }
284                        else if(tmp.getYvalue() == -2) {
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         * Return a set of Matches identified in Dynamic programming matrix. A match
332         * is a pair of subsequences whose score is higher than the preset
333         * scoreThreshold
334         **/
335        public List<Match> getMatches() {
336                ArrayList<Match> matchList = new ArrayList();
337                int fA = 0, fB = 0;
338                // skip the first row and column, find the next maxScore after
339                // prevmaxScore
340                for (int i = 1; i <= length1; i++) {
341                        for (int j = 1; j <= length2; j++) {
342                                if (matrix[i][j].getScore() > scoreThreshold
343                                                && matrix[i][j].getScore() > matrix[i - 1][j - 1].getScore()
344                                                && matrix[i][j].getScore() > matrix[i - 1][j].getScore()
345                                                && matrix[i][j].getScore() > matrix[i][j - 1].getScore()) {
346                                        if (i == length1 || j == length2
347                                                        || matrix[i][j].getScore() > matrix[i + 1][j + 1].getScore()) {
348                                                // should be lesser than prev maxScore
349                                                fA = i;
350                                                fB = j;
351                                                int[] f = traceback(fA, fB); // sets the x, y to
352                                                                                                                // startAlignment
353                                                                                                                // coordinates
354                                                System.out.println(f[0] + " " + i + " " + f[1] + " "
355                                                                + j + " " + matrix[i][j].getScore());
356                                                // TODO Add matches to matchList
357                                        }
358                                }
359                        }
360                }
361                return matchList;
362        }
363       
364
365        public List<NumberSequence> getAlignment() {
366                return alignment;
367        }
368
369        public void setAlignment(List<NumberSequence> alignment) {
370                this.alignment = alignment;
371        }
372
373}
Note: See TracBrowser for help on using the repository browser.