Context Navigation

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

Last change on this file since 1572 was 1572, checked in by rkrimmel, 10 years ago
Added parts of the PAL library, implemented UPGMA Algoritm for Feng Doolittle guide tree
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.matrix.SubstitutionMatrix;
7
8	public class SmithWaterman {
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: