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

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

Adding needleman wunsch
File size: 7.0 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 NeedlemanWunsch 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        private ArrayList<NumberSequence> alignment;
34
35        /**
36         * Substitution matrix to calculate scores
37         */
38        private SubstitutionMatrix submat;
39
40        public NeedlemanWunsch(int[] input1, int[] input2, SubstitutionMatrix submat,
41                        float threshold) {
42                this.input1 = input1;
43                this.input2 = input2;
44                length1 = input1.length;
45                length2 = input2.length;
46                this.submat = submat;
47
48                // System.out.println("Starting SmithWaterman algorithm with a "
49                // + submat.getClass() + " Substitution Matrix: " +
50                // submat.getClass().getCanonicalName());
51
52                matrix = new MatrixEntry[length1 + 1][length2 + 1];
53                alignment = new ArrayList<NumberSequence>();
54
55                for (int i = 0; i < length1+1; i++) {
56                        for (int j = 0; j < length2+1; j++) {
57                                matrix[i][j] = new MatrixEntry();
58                        }
59                }
60
61                buildMatrix();
62                traceback();
63        }
64
65        /**
66         * Compute the similarity score of substitution The position of the first
67         * character is 1. A position of 0 represents a gap.
68         *
69         * @param i
70         *            Position of the character in str1
71         * @param j
72         *            Position of the character in str2
73         * @return Cost of substitution of the character in str1 by the one in str2
74         */
75        private double similarity(int i, int j) {
76                return submat.getScore(input1[i - 1], input2[j - 1]);
77        }
78
79        /**
80         * Build the score matrix using dynamic programming.
81         */
82        private void buildMatrix() {
83                if (submat.getGapPenalty() >= 0) {
84                        throw new Error("Indel score must be negative");
85                }
86
87                // it's a gap
88                matrix[0][0].setScore(0);
89                matrix[0][0].setPrevious(null); // starting point
90
91                // the first column
92                for (int j = 1; j <= length2; j++) {
93                        matrix[0][j].setScore(j*submat.getGapPenalty());
94                        matrix[0][j].setPrevious(matrix[0][j - 1]);
95                        matrix[0][j].setYvalue(input2[j-1]);
96                        matrix[0][j].setXvalue(Constants.GAP_SYMBOL);
97                }
98                // the first row
99
100                for (int j = 1; j <= length1; j++) {
101                        matrix[j][0].setScore(j*submat.getGapPenalty());
102                        matrix[j][0].setPrevious(matrix[j - 1][0]);
103                        matrix[j][0].setXvalue(input1[j-1]);
104                        matrix[j][0].setYvalue(Constants.GAP_SYMBOL);
105                }
106
107                for (int i = 1; i <= length1; i++) {
108
109                        for (int j = 1; j <= length2; j++) {
110                                double diagScore = matrix[i - 1][j - 1].getScore()
111                                                + similarity(i, j);
112                                double upScore = matrix[i][j - 1].getScore()
113                                                + submat.getGapPenalty();
114                                double leftScore = matrix[i - 1][j].getScore()
115                                                + submat.getGapPenalty();
116
117                                matrix[i][j].setScore(Math.max(diagScore,
118                                                Math.max(upScore, Math.max(leftScore, 0))));
119
120                                // find the directions that give the maximum scores.
121                                // TODO: Multiple directions are ignored, we choose the first
122                                // maximum score
123                                // True if we had a match
124                                if (diagScore == matrix[i][j].getScore()) {
125                                        matrix[i][j].setPrevious(matrix[i - 1][j - 1]);
126                                        matrix[i][j].setXvalue(input1[i - 1]);
127                                        matrix[i][j].setYvalue(input2[j - 1]);
128                                }
129                                // true if we took an event from sequence x and not from y
130                                if (leftScore == matrix[i][j].getScore()) {
131                                        matrix[i][j].setXvalue(input1[i - 1]);
132                                        matrix[i][j].setYvalue(Constants.GAP_SYMBOL);
133                                        matrix[i][j].setPrevious(matrix[i - 1][j]);
134                                }
135                                // true if we took an event from sequence y and not from x
136                                if (upScore == matrix[i][j].getScore()) {
137                                        matrix[i][j].setXvalue(Constants.GAP_SYMBOL);
138                                        matrix[i][j].setYvalue(input2[j - 1]);
139                                        matrix[i][j].setPrevious(matrix[i][j - 1]);
140                                }
141                        }
142                }
143        }
144
145        /**
146         * Get the maximum value in the score matrix.
147         */
148        public double getMaxScore() {
149                double maxScore = 0;
150
151                // skip the first row and column
152                for (int i = 1; i <= length1; i++) {
153                        for (int j = 1; j <= length2; j++) {
154                                if (matrix[i][j].getScore() > maxScore) {
155                                        maxScore = matrix[i][j].getScore();
156                                }
157                        }
158                }
159
160                return maxScore;
161        }
162
163        /*
164         * (non-Javadoc)
165         *
166         * @see
167         * de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.AlignmentAlgorithm
168         * #getAlignmentScore()
169         */
170        @Override
171        public double getAlignmentScore() {
172                return getMaxScore();
173        }
174
175        public void traceback() {
176                MatrixEntry tmp = matrix[length1][length2];
177                int aligned1[] = new int[length1 + length2 + 2];
178                int aligned2[] = new int[length1 + length2 + 2];
179                int count = 0;
180                do {
181                        if (length1 + length2 + 2 == count) {
182                                Console.traceln(Level.WARNING,
183                                                "Traceback longer than both sequences summed up!");
184                                break;
185                        }
186                        aligned1[count] = tmp.getXvalue();
187                        aligned2[count] = tmp.getYvalue();
188
189                        tmp = tmp.getPrevious();
190                        count++;
191
192                } while (tmp != null);
193                count--;
194                // reverse order of the alignment
195                int reversed1[] = new int[count];
196                int reversed2[] = new int[count];
197
198                for (int i = count; i > 0; i--) {
199                        reversed1[reversed1.length - i] = aligned1[i];
200                        reversed2[reversed2.length - i] = aligned2[i];
201                }
202
203                NumberSequence ns1 = new NumberSequence(reversed1.length);
204                NumberSequence ns2 = new NumberSequence(reversed2.length);
205                ns1.setSequence(reversed1);
206                ns2.setSequence(reversed2);
207
208                alignment.add(ns1);
209                alignment.add(ns2);
210        }
211
212        public void printAlignment() {
213                MatrixEntry tmp = matrix[length1][length2];
214                String aligned1 = "";
215                String aligned2 = "";
216                int count = 0;
217                do {
218                        String append1 = "";
219                        String append2 = "";
220
221                        if (tmp.getXvalue() == Constants.GAP_SYMBOL) {
222                                append1 = "  ___";
223                        } else {
224                                append1 = String.format("%5d", tmp.getXvalue());
225                        }
226
227                        if (tmp.getYvalue() == Constants.GAP_SYMBOL) {
228                                append2 = "  ___";
229                        } else {
230                                append2 = String.format("%5d", tmp.getYvalue());
231                        }
232                        if (count != 0) {
233                                aligned1 = append1 + aligned1;
234                                aligned2 = append2 + aligned2;
235                        }
236
237                        tmp = tmp.getPrevious();
238                        count++;
239
240                } while (tmp != null);
241                System.out.println(aligned1);
242                System.out.println(aligned2);
243        }
244
245        /**
246         * print the dynmaic programming matrix
247         */
248        public void printDPMatrix() {
249                System.out.print("          ");
250                for (int i = 1; i <= length1; i++)
251                        System.out.format("%5d", input1[i - 1]);
252                System.out.println();
253                for (int j = 0; j <= length2; j++) {
254                        if (j > 0)
255                                System.out.format("%5d ", input2[j - 1]);
256                        else {
257                                System.out.print("      ");
258                        }
259                        for (int i = 0; i <= length1; i++) {
260                                        System.out.format("%4.1f ", matrix[i][j].getScore());
261                        }
262                        System.out.println();
263                }
264        }
265
266        /*
267         * (non-Javadoc)
268         *
269         * @see
270         * de.ugoe.cs.autoquest.tasktrees.alignment.algorithms.AlignmentAlgorithm
271         * #getAlignment()
272         */
273        @Override
274        public ArrayList<NumberSequence> getAlignment() {
275                return alignment;
276        }
277
278        public void setAlignment(ArrayList<NumberSequence> alignment) {
279                this.alignment = alignment;
280        }
281
282
283
284}
Note: See TracBrowser for help on using the repository browser.