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

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

Startet implementing signigicant match detection
File size: 9.8 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.logging.Level;
7
8import de.ugoe.cs.autoquest.tasktrees.alignment.matrix.SubstitutionMatrix;
9import de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.Constants;
10import de.ugoe.cs.util.console.Console;
11
12public class SmithWatermanRepeated implements AlignmentAlgorithm {
13
14        /**
15         * The first input
16         */
17        private int[] input1;
18
19        /**
20         * The second input String
21         */
22        private int[] input2;
23
24        /**
25         * The lengths of the input
26         */
27        private int length1, length2;
28
29        /**
30         * The score matrix. The true scores should be divided by the normalization
31         * factor.
32         */
33        private MatrixEntry[][] matrix;
34
35
36        private ArrayList<NumberSequence> alignment;
37       
38        private float scoreThreshold;
39
40        /**
41         * Substitution matrix to calculate scores
42         */
43        private SubstitutionMatrix submat;
44
45        public SmithWatermanRepeated() {
46       
47        }
48
49        /**
50         * Compute the similarity score of substitution The position of the first
51         * character is 1. A position of 0 represents a gap.
52         *
53         * @param i
54         *            Position of the character in str1
55         * @param j
56         *            Position of the character in str2
57         * @return Cost of substitution of the character in str1 by the one in str2
58         */
59        private double similarity(int i, int j) {
60                return submat.getScore(input1[i - 1], input2[j - 1]);
61        }
62
63        /**
64         * Build the score matrix using dynamic programming.
65         */
66        private void buildMatrix() {
67                if (submat.getGapPenalty() >= 0) {
68                        throw new Error("Indel score must be negative");
69                }
70               
71                // it's a gap
72                matrix[0][0].setScore(0);
73                matrix[0][0].setPrevious(null); // starting point
74                matrix[0][0].setXvalue(Constants.UNMATCHED_SYMBOL);
75                matrix[0][0].setYvalue(Constants.UNMATCHED_SYMBOL);
76
77                // the first column
78                for (int j = 1; j < length2; j++) {
79                        matrix[0][j].setScore(0);
80                        //We don't need to go back to [0][0] if we reached matrix[0][x], so just end here
81                        //matrix[0][j].setPrevious(matrix[0][j-1]);
82                        matrix[0][j].setPrevious(null);
83                }
84               
85               
86               
87                for (int i = 1; i < length1 + 2; i++) {
88               
89                        // Formula for first row:
90                        // F(i,0) = max { F(i-1,0), F(i-1,j)-T  j=1,...,m
91                       
92                        double firstRowLeftScore = matrix[i-1][0].getScore();
93                        //for sequences of length 1
94                        double tempMax;
95                        int maxRowIndex;
96                        if(length2 == 1) {
97                                tempMax = matrix[i-1][0].getScore();
98                                maxRowIndex = 0;
99                        } else {
100                                tempMax = matrix[i-1][1].getScore();
101                                maxRowIndex = 1;
102                                //position of the maximal score of the previous row
103                               
104                                for(int j = 2; j < length2;j++) {
105                                        if(matrix[i-1][j].getScore() > tempMax) {
106                                                tempMax = matrix[i-1][j].getScore();
107                                                maxRowIndex = j;
108                                        }
109                                }
110
111                        }
112                                       
113                        tempMax -= scoreThreshold;
114                        matrix[i][0].setScore(Math.max(firstRowLeftScore, tempMax));
115                        if(tempMax == matrix[i][0].getScore()){
116                                matrix[i][0].setPrevious(matrix[i-1][maxRowIndex]);
117                        }
118                       
119                        if(firstRowLeftScore == matrix[i][0].getScore()) {
120                                matrix[i][0].setPrevious(matrix[i-1][0]);
121                        }
122                       
123                        //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
124                        if(i<length1+1)
125                        {
126                                matrix[i][0].setXvalue(input1[i-1]);
127                                matrix[i][0].setYvalue(Constants.UNMATCHED_SYMBOL);
128                        }
129                        else {
130                        //End after we calculated final score
131                                return;
132                        }
133                       
134                       
135                        for (int j = 1; j < length2; j++) {
136                                double diagScore = matrix[i - 1][j - 1].getScore() + similarity(i, j);
137                                double upScore = matrix[i][j - 1].getScore() + submat.getGapPenalty();
138                                double leftScore = matrix[i - 1][j].getScore() + submat.getGapPenalty();
139
140                                matrix[i][j].setScore(Math.max(diagScore,Math.max(upScore, Math.max(leftScore,matrix[i][0].getScore()))));
141
142                                // find the directions that give the maximum scores.
143                                // TODO: Multiple directions are ignored, we choose the first maximum score
144                                //True if we had a match
145                                if (diagScore == matrix[i][j].getScore()) {
146                                        matrix[i][j].setPrevious(matrix[i-1][j-1]);
147                                        matrix[i][j].setXvalue(input1[i-1]);
148                                        matrix[i][j].setYvalue(input2[j-1]);
149                                }
150                                //true if we took an event from sequence x and not from y
151                                if (leftScore == matrix[i][j].getScore()) {
152                                        matrix[i][j].setXvalue(input1[i-1]);
153                                        matrix[i][j].setYvalue(Constants.GAP_SYMBOL);
154                                        matrix[i][j].setPrevious(matrix[i-1][j]);
155                                }
156                                //true if we took an event from sequence y and not from x
157                                if (upScore == matrix[i][j].getScore()) {
158                                        matrix[i][j].setXvalue(Constants.GAP_SYMBOL);
159                                        matrix[i][j].setYvalue(input2[j-1]);
160                                        matrix[i][j].setPrevious(matrix[i][j-1]);
161                                }
162                                //true if we ended a matching region
163                                if (matrix[i][0].getScore() == matrix[i][j].getScore()) {
164                                        matrix[i][j].setPrevious(matrix[i][0]);
165                                        matrix[i][j].setXvalue(input1[i-1]);
166                                        matrix[i][j].setYvalue(Constants.UNMATCHED_SYMBOL);
167                                }
168                        }
169                       
170                        //Set the complete score cell
171                       
172                }
173        }
174
175        /**
176         * Get the maximum value in the score matrix.
177         */
178        public double getMaxScore() {
179                double maxScore = 0;
180
181                // skip the first row and column
182                for (int i = 1; i <= length1; i++) {
183                        for (int j = 1; j < length2; j++) {
184                                if (matrix[i][j].getScore() > maxScore) {
185                                        maxScore = matrix[i][j].getScore();
186                                }
187                        }
188                }
189
190                return maxScore;
191        }
192
193        /* (non-Javadoc)
194         * @see de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.AlignmentAlgorithm#getAlignmentScore()
195         */
196        @Override
197        public double getAlignmentScore() {
198                return matrix[length1+1][0].getScore();
199        }
200
201        public void traceback() {
202                MatrixEntry tmp = matrix[length1+1][0].getPrevious();
203                LinkedList<Integer> aligned1 = new LinkedList<Integer>();
204                LinkedList<Integer> aligned2 = new LinkedList<Integer>();
205                while (tmp.getPrevious() != null) {
206                       
207                        aligned1.add(new Integer(tmp.getXvalue()));
208                        aligned2.add(new Integer(tmp.getYvalue()));
209
210                        tmp = tmp.getPrevious();
211                }
212               
213                // reverse order of the alignment
214                int reversed1[] = new int[aligned1.size()];
215                int reversed2[] = new int[aligned2.size()];
216
217                int count = 0;
218                for (Iterator<Integer> it = aligned1.iterator(); it.hasNext();) {
219                        count++;
220                        reversed1[reversed1.length - count] = it.next();
221                       
222                }
223                count = 0;
224                for (Iterator<Integer> it = aligned2.iterator(); it.hasNext();) {
225                        count++;
226                        reversed2[reversed2.length - count] = it.next();
227                }
228
229                NumberSequence ns1 = new NumberSequence(reversed1.length);
230                NumberSequence ns2 = new NumberSequence(reversed2.length);
231                ns1.setSequence(reversed1);
232                ns2.setSequence(reversed2);
233
234                alignment.add(ns1);
235                alignment.add(ns2);
236        }
237       
238
239       
240        public void printAlignment() {
241                int[] tmp1 = alignment.get(0).getSequence();
242                int[] tmp2 = alignment.get(1).getSequence();
243                for (int i=0; i< tmp1.length;i++) {
244                        if(tmp1[i] == Constants.GAP_SYMBOL) {
245                                System.out.print("  ___");
246                        }
247                        else if(tmp1[i] == Constants.UNMATCHED_SYMBOL) {
248                                System.out.print("  ...");
249                        }
250                        else {
251                                System.out.format("%5d", tmp1[i]);
252                        }
253                       
254                }
255                System.out.println();
256                for (int i=0; i< tmp2.length;i++) {
257                        if(tmp2[i] == Constants.GAP_SYMBOL) {
258                                System.out.print("  ___");
259                        }
260                        else if(tmp2[i] == Constants.UNMATCHED_SYMBOL) {
261                                System.out.print("  ...");
262                        }
263                        else {
264                                System.out.format("%5d", tmp2[i]);
265                        }
266                       
267                }
268                System.out.println();
269               
270               
271       
272        }
273       
274        public ArrayList<ArrayList<NumberSequence>> getMatches() {
275                ArrayList<ArrayList<NumberSequence>> result = new ArrayList<ArrayList<NumberSequence>>();
276               
277                //both alignment sequences should be equally long
278                int i = 0;
279                int[] seq1 = alignment.get(0).getSequence();
280                int[] seq2 = alignment.get(1).getSequence();
281                int start = 0;
282                while (i < seq1.length){
283                        if(seq2[i] != Constants.UNMATCHED_SYMBOL) {
284                                start = i;
285                                int count = 0;
286                                while(i < seq2.length && seq2[i] != Constants.UNMATCHED_SYMBOL) {
287                                        i++;
288                                        count++;
289                                }
290                                //I am really missing memcpy here
291                                int[] tmp1 = new int[count];
292                                int[] tmp2 = new int[count];
293                                for (int j = 0; j<count;j++) {
294                                        tmp1[j] = seq1[start+j];
295                                        tmp2[j] = seq2[start+j];
296                                }
297                                NumberSequence tmpns1 = new NumberSequence(count);
298                                NumberSequence tmpns2 = new NumberSequence(count);
299                                tmpns1.setSequence(tmp1);
300                                tmpns2.setSequence(tmp2);
301                                ArrayList<NumberSequence> tmpal = new ArrayList<NumberSequence>();
302                                tmpal.add(tmpns1);
303                                tmpal.add(tmpns2);
304                                result.add(tmpal);
305                        }
306                        i++;
307                }
308               
309               
310               
311               
312               
313               
314                return result;
315               
316        }
317       
318        /**
319         * print the dynmaic programming matrix
320         */
321        public void printDPMatrix() {
322                System.out.print("          ");
323                for (int i = 1; i <= length1; i++)
324                        System.out.format("%5d", input1[i - 1]);
325                System.out.println();
326                for (int j = 0; j <= length2; j++) {
327                        if (j > 0)
328                                System.out.format("%5d ",input2[j - 1]);
329                        else{
330                                System.out.print("      ");
331                        }
332                        for (int i = 0; i <= length1 + 1; i++) {
333                                if((i<length1+1) || (i==length1+1 && j==0))     {
334                                        System.out.format("%4.1f ",matrix[i][j].getScore());
335                                }
336                               
337                        }                       
338                        System.out.println();
339                }
340        }
341
342
343        /* (non-Javadoc)
344         * @see de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.AlignmentAlgorithm#getAlignment()
345         */
346        @Override
347        public ArrayList<NumberSequence> getAlignment() {
348                return alignment;
349        }
350
351        public void setAlignment(ArrayList<NumberSequence> alignment) {
352                this.alignment = alignment;
353        }
354
355        @Override
356        public void align(int[] input1, int[] input2, SubstitutionMatrix submat,
357                        float threshold) {
358                this.input1 = input1;
359                this.input2 = input2;
360                length1 = input1.length;
361                length2 = input2.length;
362                this.submat = submat;
363
364                //System.out.println("Starting SmithWaterman algorithm with a "
365                //              + submat.getClass() + " Substitution Matrix: " + submat.getClass().getCanonicalName());
366                this.scoreThreshold = threshold;
367               
368                matrix = new MatrixEntry[length1+2][length2+1];
369                alignment = new ArrayList<NumberSequence>();
370               
371                for (int i = 0; i <= length1+1; i++) {
372                        for(int j = 0; j< length2; j++) {
373                                matrix[i][j] = new MatrixEntry();
374                        }
375                }
376       
377                //Console.traceln(Level.INFO,"Generating DP Matrix");
378                buildMatrix();
379                //Console.traceln(Level.INFO,"Doing traceback");
380                traceback();
381        }
382
383}
Note: See TracBrowser for help on using the repository browser.