Context Navigation

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

Last change on this file since 1558 was 1558, checked in by rkrimmel, 10 years ago
new Smith waterman variant
File size: 8.3 KB

Line
1	package de.ugoe.cs.autoquest.tasktrees.alignment.algorithms;
2
3	import java.util.ArrayList;
4	import java.util.List;
5
6	import de.ugoe.cs.autoquest.tasktrees.alignment.substitution.SubstitutionMatrix;
7
8	public class SmithWaterman 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 double[][] score;
30
31	/**
32	* Score threshold
33	*/
34	private int scoreThreshold;;
35
36	/**
37	* The similarity function constants.
38	*/
39	// static final int MATCH_SCORE = 10;
40	// static final int MISMATCH_SCORE = -8;
41	// static final int INDEL_SCORE = -9;
42
43	/**
44	* Constants of directions. Multiple directions are stored by bits. The zero
45	* direction is the starting point.
46	*/
47	static final int DR_LEFT = 1; // 0001
48	static final int DR_UP = 2; // 0010
49	static final int DR_DIAG = 4; // 0100
50	static final int DR_ZERO = 8; // 1000
51
52	/**
53	* The directions pointing to the cells that give the maximum score at the
54	* current cell. The first index is the column index. The second index is
55	* the row index.
56	*/
57	private int[][] prevCells;
58
59	/**
60	* Substitution matrix to calculate scores
61	*/
62	private SubstitutionMatrix submat;
63
64	public SmithWaterman(int[] input1, int[] input2, SubstitutionMatrix submat) {
65	this.input1 = input1;
66	this.input2 = input2;
67	length1 = input1.length;
68	length2 = input2.length;
69	this.submat = submat;
70
71	//System.out.println("Starting SmithWaterman algorithm with a "
72	// + submat.getClass() + " Substitution Matrix: " + submat.getClass().getCanonicalName());
73	scoreThreshold = 20;
74	score = new double[length1 + 1][length2 + 1];
75	prevCells = new int[length1 + 1][length2 + 1];
76
77	buildMatrix();
78	}
79
80	/**
81	* Compute the similarity score of substitution The position of the first
82	* character is 1. A position of 0 represents a gap.
83	*
84	* @param i
85	* Position of the character in str1
86	* @param j
87	* Position of the character in str2
88	* @return Cost of substitution of the character in str1 by the one in str2
89	*/
90	private double similarity(int i, int j) {
91	if (i == 0 \|\| j == 0) {
92	// it's a gap
93	return submat.getGapPenalty();
94	}
95	// System.out.println("Diag letters: " + input1[i-1] + " " +
96	// input2[j-1]);
97	// return (input1[i - 1] == input2[j - 1]) ? MATCH_SCORE :
98	// MISMATCH_SCORE;
99	return submat.getDistance(input1[i - 1], input2[j - 1]);
100	}
101
102	/**
103	* Build the score matrix using dynamic programming. Note: The indel scores
104	* must be negative. Otherwise, the part handling the first row and column
105	* has to be modified.
106	*/
107	private void buildMatrix() {
108	if (submat.getGapPenalty() >= 0) {
109	throw new Error("Indel score must be negative");
110	}
111
112	int i; // length of prefix substring of str1
113	int j; // length of prefix substring of str2
114
115	// base case
116	score[0][0] = 0;
117	prevCells[0][0] = DR_ZERO; // starting point
118
119	// the first row
120	for (i = 1; i <= length1; i++) {
121	score[i][0] = 0;
122	prevCells[i][0] = DR_ZERO;
123	}
124
125	// the first column
126	for (j = 1; j <= length2; j++) {
127	score[0][j] = 0;
128	prevCells[0][j] = DR_ZERO;
129	}
130
131	// the rest of the matrix
132	for (i = 1; i <= length1; i++) {
133	for (j = 1; j <= length2; j++) {
134	double diagScore = score[i - 1][j - 1] + similarity(i, j);
135	double upScore = score[i][j - 1] + similarity(0, j);
136	double leftScore = score[i - 1][j] + similarity(i, 0);
137
138	score[i][j] = Math.max(diagScore,
139	Math.max(upScore, Math.max(leftScore, 0)));
140	prevCells[i][j] = 0;
141
142	// find the directions that give the maximum scores.
143	// the bitwise OR operator is used to record multiple
144	// directions.
145	if (diagScore == score[i][j]) {
146	prevCells[i][j] \|= DR_DIAG;
147	}
148	if (leftScore == score[i][j]) {
149	prevCells[i][j] \|= DR_LEFT;
150	}
151	if (upScore == score[i][j]) {
152	prevCells[i][j] \|= DR_UP;
153	}
154	if (0 == score[i][j]) {
155	prevCells[i][j] \|= DR_ZERO;
156	}
157	}
158	}
159	}
160
161	/**
162	* Get the maximum value in the score matrix.
163	*/
164	public double getMaxScore() {
165	double maxScore = 0;
166
167	// skip the first row and column
168	for (int i = 1; i <= length1; i++) {
169	for (int j = 1; j <= length2; j++) {
170	if (score[i][j] > maxScore) {
171	maxScore = score[i][j];
172	}
173	}
174	}
175
176	return maxScore;
177	}
178
179	/**
180	* Get the alignment score between the two input strings.
181	*/
182	public double getAlignmentScore() {
183	return getMaxScore();
184	}
185
186
187
188	/**
189	* TODO: Iterative Version!!! Output the local alignments ending in the (i,
190	* j) cell. aligned1 and aligned2 are suffixes of final aligned strings
191	* found in backtracking before calling this function. Note: the strings are
192	* replicated at each recursive call. Use buffers or stacks to improve
193	* efficiency.
194	*/
195	private void printAlignments(int i, int j, String aligned1, String aligned2) {
196	// we've reached the starting point, so print the alignments
197
198	if ((prevCells[i][j] & DR_ZERO) > 0) {
199	System.out.println(aligned1);
200	System.out.println(aligned2);
201	System.out.println();
202
203	// Note: we could check other directions for longer alignments
204	// with the same score. we don't do it here.
205	return;
206	}
207
208	// find out which directions to backtrack
209	if ((prevCells[i][j] & DR_LEFT) > 0) {
210	printAlignments(i - 1, j, input1[i - 1] + aligned1, "_ " + aligned2);
211	}
212	if ((prevCells[i][j] & DR_UP) > 0) {
213	printAlignments(i, j - 1, "_ " + aligned1, input2[j - 1] + aligned2);
214	}
215	if ((prevCells[i][j] & DR_DIAG) > 0) {
216	printAlignments(i - 1, j - 1, input1[i - 1] + " " + aligned1,
217	input2[j - 1] + " " + aligned2);
218	}
219	}
220
221	/**
222	* given the bottom right corner point trace back the top left conrner. at
223	* entry: i, j hold bottom right (end of Aligment coords) at return: hold
224	* top left (start of Alignment coords)
225	*/
226	private int[] traceback(int i, int j) {
227
228	// find out which directions to backtrack
229	while (true) {
230	if ((prevCells[i][j] & DR_LEFT) > 0) {
231	if (score[i - 1][j] > 0)
232	i--;
233	else
234	break;
235	}
236	if ((prevCells[i][j] & DR_UP) > 0) {
237	// return traceback(i, j-1);
238	if (score[i][j - 1] > 0)
239	j--;
240	else
241	break;
242	}
243	if ((prevCells[i][j] & DR_DIAG) > 0) {
244	// return traceback(i-1, j-1);
245	if (score[i - 1][j - 1] > 0) {
246	i--;
247	j--;
248	} else
249	break;
250	}
251	}
252	int[] m = { i, j };
253	return m;
254	}
255
256	/**
257	* Output the local alignments with the maximum score.
258	*/
259	public void printAlignments() {
260	// find the cell with the maximum score
261	double maxScore = getMaxScore();
262
263	/*
264	* for (int i = 0; i < matches.length; i++) {
265	* System.out.println("Match #" + i + ":" + matches.get(i)); }
266	*/
267
268	// skip the first row and column
269	for (int i = 1; i <= length1; i++) {
270	for (int j = 1; j <= length2; j++) {
271	if (score[i][j] == maxScore) {
272	printAlignments(i, j, "", "");
273	}
274	}
275	}
276	}
277
278	/**
279	* print the dynmaic programming matrix
280	*/
281	public void printDPMatrix() {
282	System.out.print(" ");
283	for (int j = 1; j <= length2; j++)
284	System.out.format("%5d", input2[j - 1]);
285	System.out.println();
286	for (int i = 0; i <= length1; i++) {
287	if (i > 0)
288	System.out.format("%5d ",input1[i - 1]);
289	else{
290	System.out.print(" ");
291	}
292	for (int j = 0; j <= length2; j++) {
293	System.out.format("%4.1f ",score[i][j]);
294	}
295	System.out.println();
296	}
297	}
298
299	/**
300	* Return a set of Matches identified in Dynamic programming matrix. A match
301	* is a pair of subsequences whose score is higher than the preset
302	* scoreThreshold
303	**/
304	public List<Match> getMatches() {
305	ArrayList<Match> matchList = new ArrayList();
306	int fA = 0, fB = 0;
307	// skip the first row and column, find the next maxScore after
308	// prevmaxScore
309	for (int i = 1; i <= length1; i++) {
310	for (int j = 1; j <= length2; j++) {
311	if (score[i][j] > scoreThreshold
312	&& score[i][j] > score[i - 1][j - 1]
313	&& score[i][j] > score[i - 1][j]
314	&& score[i][j] > score[i][j - 1]) {
315	if (i == length1 \|\| j == length2
316	\|\| score[i][j] > score[i + 1][j + 1]) {
317	// should be lesser than prev maxScore
318	fA = i;
319	fB = j;
320	int[] f = traceback(fA, fB); // sets the x, y to
321	// startAlignment
322	// coordinates
323	System.out.println(f[0] + " " + i + " " + f[1] + " "
324	+ j + " " + score[i][j]);
325	// TODO Add matches to matchList
326	}
327	}
328	}
329	}
330	return matchList;
331	}
332
333	}

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

Download in other formats: