Context Navigation

SmithWatermanRepeated.java @ 1760

Last change on this file since 1760 was 1760, checked in by rkrimmel, 10 years ago
Memory improvement and bugs
File size: 11.1 KB

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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: branches/autoquest-core-tasktrees-alignment/src/main/java/de/ugoe/cs/autoquest/tasktrees/alignment/algorithms/SmithWatermanRepeated.java @ 1760

Download in other formats: