Context Navigation

source: trunk/autoquest-plugin-guitar/src/main/java/de/ugoe/cs/autoquest/plugin/alignment/SmithWaterman.java @ 1303

Last change on this file since 1303 was 1303, checked in by rkrimmel, 11 years ago
Adding Match class, Smith Waterman algorithm and Substitution Matrix class
File size: 8.0 KB

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

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

Download in other formats: