Context Navigation

nwdasimple2.cpp

Visit:

Last change on this file was 10390, checked in by aboeckma, 11 years ago
added muscle sourcles amd makefile
File size: 14.0 KB

Line
1	#include "muscle.h"
2	#include "pwpath.h"
3	#include "profile.h"
4
5	#if DOUBLE_AFFINE
6
7	#define TRACE 0
8
9	extern bool g_bKeepSimpleDP;
10	extern SCORE *g_DPM;
11	extern SCORE *g_DPD;
12	extern SCORE *g_DPE;
13	extern SCORE *g_DPI;
14	extern SCORE *g_DPJ;
15	extern char *g_TBM;
16	extern char *g_TBD;
17	extern char *g_TBE;
18	extern char *g_TBI;
19	extern char *g_TBJ;
20
21	static char XlatEdgeType(char c)
22	{
23	if ('E' == c)
24	return 'D';
25	if ('J' == c)
26	return 'I';
27	return c;
28	}
29
30	static const char *LocalScoreToStr(SCORE s)
31	{
32	static char str[16];
33	if (s < -100000)
34	return " *";
35	sprintf(str, "%6.1f", s);
36	return str;
37	}
38
39	static void ListDP(const SCORE DPM_, const ProfPos PA, const ProfPos *PB,
40	unsigned uPrefixCountA, unsigned uPrefixCountB)
41	{
42	Log(" ");
43	for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
44	{
45	char c = ' ';
46	if (uPrefixLengthB > 0)
47	c = ConsensusChar(PB[uPrefixLengthB - 1]);
48	Log(" %4u:%c", uPrefixLengthB, c);
49	}
50	Log("\n");
51	for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
52	{
53	char c = ' ';
54	if (uPrefixLengthA > 0)
55	c = ConsensusChar(PA[uPrefixLengthA - 1]);
56	Log("%4u:%c ", uPrefixLengthA, c);
57	for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
58	Log(" %s", LocalScoreToStr(DPM(uPrefixLengthA, uPrefixLengthB)));
59	Log("\n");
60	}
61	}
62
63	static void ListTB(const char TBM_, const ProfPos PA, const ProfPos *PB,
64	unsigned uPrefixCountA, unsigned uPrefixCountB)
65	{
66	Log(" ");
67	for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
68	{
69	char c = ' ';
70	if (uPrefixLengthB > 0)
71	c = ConsensusChar(PB[uPrefixLengthB - 1]);
72	Log(" %4u:%c", uPrefixLengthB, c);
73	}
74	Log("\n");
75	for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
76	{
77	char c = ' ';
78	if (uPrefixLengthA > 0)
79	c = ConsensusChar(PA[uPrefixLengthA - 1]);
80	Log("%4u:%c ", uPrefixLengthA, c);
81	for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
82	Log(" %6c", TBM(uPrefixLengthA, uPrefixLengthB));
83	Log("\n");
84	}
85	}
86
87	static void ListDPM(const SCORE DPM_, const ProfPos PA, const ProfPos *PB,
88	unsigned uPrefixCountA, unsigned uPrefixCountB)
89	{
90	Log(" ");
91	for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
92	{
93	char c = ' ';
94	if (uPrefixLengthB > 0)
95	c = ConsensusChar(PB[uPrefixLengthB - 1]);
96	Log(" %4u:%c", uPrefixLengthB, c);
97	}
98	Log("\n");
99	for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
100	{
101	char c = ' ';
102	if (uPrefixLengthA > 0)
103	c = ConsensusChar(PA[uPrefixLengthA - 1]);
104	Log("%4u:%c ", uPrefixLengthA, c);
105	for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
106	{
107	SCORE x = (uPrefixLengthA + uPrefixLengthB)*g_scoreGapExtend;
108	SCORE s = DPM(uPrefixLengthA, uPrefixLengthB) - x;
109	Log(" %s", LocalScoreToStr(s));
110	}
111	Log("\n");
112	}
113	}
114
115	extern SCORE ScoreProfPos2(const ProfPos &PP, const ProfPos &PPB);
116
117	SCORE NWDASimple2(const ProfPos PA, unsigned uLengthA, const ProfPos PB,
118	unsigned uLengthB, PWPath &Path)
119	{
120	assert(uLengthB > 0 && uLengthA > 0);
121
122	const unsigned uPrefixCountA = uLengthA + 1;
123	const unsigned uPrefixCountB = uLengthB + 1;
124
125	// Allocate DP matrices
126	const size_t LM = uPrefixCountA*uPrefixCountB;
127	SCORE *DPM_ = new SCORE[LM];
128	SCORE *DPD_ = new SCORE[LM];
129	SCORE *DPE_ = new SCORE[LM];
130	SCORE *DPI_ = new SCORE[LM];
131	SCORE *DPJ_ = new SCORE[LM];
132	SCORE *DPL_ = new SCORE[LM];
133
134	char *TBM_ = new char[LM];
135	char *TBD_ = new char[LM];
136	char *TBE_ = new char[LM];
137	char *TBI_ = new char[LM];
138	char *TBJ_ = new char[LM];
139
140	memset(DPM_, 0, LM*sizeof(SCORE));
141	memset(DPD_, 0, LM*sizeof(SCORE));
142	memset(DPE_, 0, LM*sizeof(SCORE));
143	memset(DPI_, 0, LM*sizeof(SCORE));
144	memset(DPJ_, 0, LM*sizeof(SCORE));
145
146	// memset(DPL_, 0, LM*sizeof(SCORE));
147
148	memset(TBM_, '?', LM);
149	memset(TBD_, '?', LM);
150	memset(TBE_, '?', LM);
151	memset(TBI_, '?', LM);
152	memset(TBJ_, '?', LM);
153
154	DPM(0, 0) = 0;
155	DPD(0, 0) = MINUS_INFINITY;
156	DPE(0, 0) = MINUS_INFINITY;
157	DPI(0, 0) = MINUS_INFINITY;
158	DPJ(0, 0) = MINUS_INFINITY;
159
160	DPM(1, 0) = MINUS_INFINITY;
161	DPD(1, 0) = PA[0].m_scoreGapOpen;
162	DPE(1, 0) = PA[0].m_scoreGapOpen2;
163	DPI(1, 0) = MINUS_INFINITY;
164	DPJ(1, 0) = MINUS_INFINITY;
165
166	DPM(0, 1) = MINUS_INFINITY;
167	DPD(0, 1) = MINUS_INFINITY;
168	DPE(0, 1) = MINUS_INFINITY;
169	DPI(0, 1) = PB[0].m_scoreGapOpen;
170	DPJ(0, 1) = PB[0].m_scoreGapOpen2;
171
172	// Empty prefix of B is special case
173	for (unsigned uPrefixLengthA = 2; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
174	{
175	// M=LetterA+LetterB, impossible with empty prefix
176	DPM(uPrefixLengthA, 0) = MINUS_INFINITY;
177
178	// D=LetterA+GapB
179	DPD(uPrefixLengthA, 0) = DPD(uPrefixLengthA - 1, 0) + g_scoreGapExtend;
180	TBD(uPrefixLengthA, 0) = 'D';
181
182	DPE(uPrefixLengthA, 0) = DPE(uPrefixLengthA - 1, 0) + g_scoreGapExtend2;
183	TBE(uPrefixLengthA, 0) = 'E';
184
185	// I=GapA+LetterB, impossible with empty prefix
186	DPI(uPrefixLengthA, 0) = MINUS_INFINITY;
187	DPJ(uPrefixLengthA, 0) = MINUS_INFINITY;
188	}
189
190	// Empty prefix of A is special case
191	for (unsigned uPrefixLengthB = 2; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
192	{
193	// M=LetterA+LetterB, impossible with empty prefix
194	DPM(0, uPrefixLengthB) = MINUS_INFINITY;
195
196	// D=LetterA+GapB, impossible with empty prefix
197	DPD(0, uPrefixLengthB) = MINUS_INFINITY;
198	DPE(0, uPrefixLengthB) = MINUS_INFINITY;
199
200	// I=GapA+LetterB
201	DPI(0, uPrefixLengthB) = DPI(0, uPrefixLengthB - 1) + g_scoreGapExtend;
202	TBI(0, uPrefixLengthB) = 'I';
203
204	DPJ(0, uPrefixLengthB) = DPJ(0, uPrefixLengthB - 1) + g_scoreGapExtend2;
205	TBJ(0, uPrefixLengthB) = 'J';
206	}
207
208	// ============
209	// Main DP loop
210	// ============
211	for (unsigned uPrefixLengthB = 1; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
212	{
213	const ProfPos &PPB = PB[uPrefixLengthB - 1];
214	SCORE scoreGapCloseB;
215	if (uPrefixLengthB == 1)
216	scoreGapCloseB = MINUS_INFINITY;
217	else
218	scoreGapCloseB = PB[uPrefixLengthB-2].m_scoreGapClose;
219
220	SCORE scoreGapClose2B;
221	if (uPrefixLengthB == 1)
222	scoreGapClose2B = MINUS_INFINITY;
223	else
224	scoreGapClose2B = PB[uPrefixLengthB-2].m_scoreGapClose2;
225
226	for (unsigned uPrefixLengthA = 1; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
227	{
228	const ProfPos &PPA = PA[uPrefixLengthA - 1];
229
230	{
231	// Match M=LetterA+LetterB
232	SCORE scoreLL = ScoreProfPos2(PPA, PPB);
233	DPL(uPrefixLengthA, uPrefixLengthB) = scoreLL;
234
235	SCORE scoreGapCloseA;
236	if (uPrefixLengthA == 1)
237	scoreGapCloseA = MINUS_INFINITY;
238	else
239	scoreGapCloseA = PA[uPrefixLengthA-2].m_scoreGapClose;
240
241	SCORE scoreGapClose2A;
242	if (uPrefixLengthA == 1)
243	scoreGapClose2A = MINUS_INFINITY;
244	else
245	scoreGapClose2A = PA[uPrefixLengthA-2].m_scoreGapClose2;
246
247	SCORE scoreMM = DPM(uPrefixLengthA-1, uPrefixLengthB-1);
248	SCORE scoreDM = DPD(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseA;
249	SCORE scoreEM = DPE(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapClose2A;
250	SCORE scoreIM = DPI(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseB;
251	SCORE scoreJM = DPJ(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapClose2B;
252	SCORE scoreBest;
253	if (scoreMM >= scoreDM && scoreMM >= scoreIM && scoreMM >= scoreEM && scoreMM >= scoreJM)
254	{
255	scoreBest = scoreMM;
256	TBM(uPrefixLengthA, uPrefixLengthB) = 'M';
257	}
258	else if (scoreDM >= scoreMM && scoreDM >= scoreIM && scoreDM >= scoreEM && scoreDM >= scoreJM)
259	{
260	scoreBest = scoreDM;
261	TBM(uPrefixLengthA, uPrefixLengthB) = 'D';
262	}
263	else if (scoreEM >= scoreMM && scoreEM >= scoreIM && scoreEM >= scoreDM && scoreEM >= scoreJM)
264	{
265	scoreBest = scoreEM;
266	TBM(uPrefixLengthA, uPrefixLengthB) = 'E';
267	}
268	else if (scoreIM >= scoreMM && scoreIM >= scoreDM && scoreIM >= scoreEM && scoreIM >= scoreJM)
269	{
270	scoreBest = scoreIM;
271	TBM(uPrefixLengthA, uPrefixLengthB) = 'I';
272	}
273	else if (scoreJM >= scoreMM && scoreJM >= scoreDM && scoreJM >= scoreEM && scoreJM >= scoreIM)
274	{
275	scoreBest = scoreJM;
276	TBM(uPrefixLengthA, uPrefixLengthB) = 'J';
277	}
278	else
279	Quit("Max failed (M)");
280
281	DPM(uPrefixLengthA, uPrefixLengthB) = scoreBest + scoreLL;
282	}
283
284	{
285	// Delete D=LetterA+GapB
286	SCORE scoreMD = DPM(uPrefixLengthA-1, uPrefixLengthB) +
287	PA[uPrefixLengthA-1].m_scoreGapOpen;
288	SCORE scoreDD = DPD(uPrefixLengthA-1, uPrefixLengthB) +
289	g_scoreGapExtend;
290
291	SCORE scoreBest;
292	if (scoreMD >= scoreDD)
293	{
294	scoreBest = scoreMD;
295	TBD(uPrefixLengthA, uPrefixLengthB) = 'M';
296	}
297	else
298	{
299	assert(scoreDD >= scoreMD);
300	scoreBest = scoreDD;
301	TBD(uPrefixLengthA, uPrefixLengthB) = 'D';
302	}
303	DPD(uPrefixLengthA, uPrefixLengthB) = scoreBest;
304	}
305
306	{
307	// Delete E=LetterA+GapB
308	SCORE scoreME = DPM(uPrefixLengthA-1, uPrefixLengthB) +
309	PA[uPrefixLengthA-1].m_scoreGapOpen2;
310	SCORE scoreEE = DPE(uPrefixLengthA-1, uPrefixLengthB) +
311	g_scoreGapExtend2;
312
313	SCORE scoreBest;
314	if (scoreME >= scoreEE)
315	{
316	scoreBest = scoreME;
317	TBE(uPrefixLengthA, uPrefixLengthB) = 'M';
318	}
319	else
320	{
321	assert(scoreEE >= scoreME);
322	scoreBest = scoreEE;
323	TBE(uPrefixLengthA, uPrefixLengthB) = 'E';
324	}
325	DPE(uPrefixLengthA, uPrefixLengthB) = scoreBest;
326	}
327
328	// Insert I=GapA+LetterB
329	{
330	SCORE scoreMI = DPM(uPrefixLengthA, uPrefixLengthB-1) +
331	PB[uPrefixLengthB-1].m_scoreGapOpen;
332	SCORE scoreII = DPI(uPrefixLengthA, uPrefixLengthB-1) +
333	g_scoreGapExtend;
334
335	SCORE scoreBest;
336	if (scoreMI >= scoreII)
337	{
338	scoreBest = scoreMI;
339	TBI(uPrefixLengthA, uPrefixLengthB) = 'M';
340	}
341	else
342	{
343	assert(scoreII > scoreMI);
344	scoreBest = scoreII;
345	TBI(uPrefixLengthA, uPrefixLengthB) = 'I';
346	}
347	DPI(uPrefixLengthA, uPrefixLengthB) = scoreBest;
348	}
349
350	// Insert J=GapA+LetterB
351	{
352	SCORE scoreMJ = DPM(uPrefixLengthA, uPrefixLengthB-1) +
353	PB[uPrefixLengthB-1].m_scoreGapOpen2;
354	SCORE scoreJJ = DPJ(uPrefixLengthA, uPrefixLengthB-1) +
355	g_scoreGapExtend2;
356
357	SCORE scoreBest;
358	if (scoreMJ > scoreJJ)
359	{
360	scoreBest = scoreMJ;
361	TBJ(uPrefixLengthA, uPrefixLengthB) = 'M';
362	}
363	else
364	{
365	assert(scoreJJ >= scoreMJ);
366	scoreBest = scoreJJ;
367	TBJ(uPrefixLengthA, uPrefixLengthB) = 'J';
368	}
369	DPJ(uPrefixLengthA, uPrefixLengthB) = scoreBest;
370	}
371	}
372	}
373
374	// Special case: close gaps at end of alignment
375	DPD(uLengthA, uLengthB) += PA[uLengthA-1].m_scoreGapClose;
376	DPE(uLengthA, uLengthB) += PA[uLengthA-1].m_scoreGapClose2;
377
378	DPI(uLengthA, uLengthB) += PB[uLengthB-1].m_scoreGapClose;
379	DPJ(uLengthA, uLengthB) += PB[uLengthB-1].m_scoreGapClose2;
380
381	#if TRACE
382	Log("DPL:\n");
383	ListDP(DPL_, PA, PB, uPrefixCountA, uPrefixCountB);
384	Log("DPM:\n");
385	ListDP(DPM_, PA, PB, uPrefixCountA, uPrefixCountB);
386	Log("DPD:\n");
387	ListDP(DPD_, PA, PB, uPrefixCountA, uPrefixCountB);
388	Log("DPE:\n");
389	ListDP(DPE_, PA, PB, uPrefixCountA, uPrefixCountB);
390	Log("DPI:\n");
391	ListDP(DPI_, PA, PB, uPrefixCountA, uPrefixCountB);
392	Log("DPJ:\n");
393	ListDP(DPJ_, PA, PB, uPrefixCountA, uPrefixCountB);
394	Log("TBM:\n");
395	ListTB(TBM_, PA, PB, uPrefixCountA, uPrefixCountB);
396	Log("TBD:\n");
397	ListTB(TBD_, PA, PB, uPrefixCountA, uPrefixCountB);
398	Log("TBE:\n");
399	ListTB(TBE_, PA, PB, uPrefixCountA, uPrefixCountB);
400	Log("TBI:\n");
401	ListTB(TBI_, PA, PB, uPrefixCountA, uPrefixCountB);
402	Log("TBJ:\n");
403	ListTB(TBJ_, PA, PB, uPrefixCountA, uPrefixCountB);
404	#endif
405
406	// ==========
407	// Trace-back
408	// ==========
409
410	Path.Clear();
411
412	// Find last edge
413	char cEdgeType = '?';
414	SCORE BestScore = MINUS_INFINITY;
415	SCORE M = DPM(uLengthA, uLengthB);
416	SCORE D = DPD(uLengthA, uLengthB);
417	SCORE E = DPE(uLengthA, uLengthB);
418	SCORE I = DPI(uLengthA, uLengthB);
419	SCORE J = DPJ(uLengthA, uLengthB);
420
421	if (M >= D && M >= E && M >= I && M >= J)
422	{
423	cEdgeType = 'M';
424	BestScore = M;
425	}
426	else if (D >= M && D >= E && D >= I && D >= J)
427	{
428	cEdgeType = 'D';
429	BestScore = D;
430	}
431	else if (E >= M && E >= D && E >= I && E >= J)
432	{
433	cEdgeType = 'E';
434	BestScore = E;
435	}
436	else if (I >= M && I >= D && I >= E && I >= J)
437	{
438	cEdgeType = 'I';
439	BestScore = I;
440	}
441	else if (J >= M && J >= D && J >= E && J >= I)
442	{
443	cEdgeType = 'J';
444	BestScore = J;
445	}
446	else
447	Quit("Bad max");
448
449	unsigned PLA = uLengthA;
450	unsigned PLB = uLengthB;
451	unsigned ECount = 0;
452	unsigned JCount = 0;
453	for (;;)
454	{
455	#if TRACE
456	Log("TraceBack: %c%u.%u\n", cEdgeType, PLA, PLB);
457	#endif
458	PWEdge Edge;
459	Edge.cType = XlatEdgeType(cEdgeType);
460	Edge.uPrefixLengthA = PLA;
461	Edge.uPrefixLengthB = PLB;
462	Path.PrependEdge(Edge);
463
464	switch (cEdgeType)
465	{
466	case 'M':
467	assert(PLA > 0);
468	assert(PLB > 0);
469	cEdgeType = TBM(PLA, PLB);
470	--PLA;
471	--PLB;
472	break;
473
474	case 'D':
475	assert(PLA > 0);
476	cEdgeType = TBD(PLA, PLB);
477	--PLA;
478	break;
479
480	case 'E':
481	++ECount;
482	assert(PLA > 0);
483	cEdgeType = TBE(PLA, PLB);
484	--PLA;
485	break;
486
487	case 'I':
488	assert(PLB > 0);
489	cEdgeType = TBI(PLA, PLB);
490	--PLB;
491	break;
492
493	case 'J':
494	++JCount;
495	assert(PLB > 0);
496	cEdgeType = TBJ(PLA, PLB);
497	--PLB;
498	break;
499
500	default:
501	Quit("Invalid edge %c", cEdgeType);
502	}
503	if (0 == PLA && 0 == PLB)
504	break;
505	}
506	//if (ECount > 0 \|\| JCount > 0)
507	// fprintf(stderr, "E=%d J=%d\n", ECount, JCount);
508	Path.Validate();
509	if (Path.GetMatchCount() + Path.GetDeleteCount() != uLengthA)
510	Quit("Path count A");
511	if (Path.GetMatchCount() + Path.GetInsertCount() != uLengthB)
512	Quit("Path count B");
513
514	if (g_bKeepSimpleDP)
515	{
516	g_DPM = DPM_;
517	g_DPD = DPD_;
518	g_DPE = DPE_;
519	g_DPI = DPI_;
520	g_DPJ = DPJ_;
521
522	g_TBM = TBM_;
523	g_TBD = TBD_;
524	g_TBE = TBE_;
525	g_TBI = TBI_;
526	g_TBJ = TBJ_;
527	}
528	else
529	{
530	delete[] DPM_;
531	delete[] DPD_;
532	delete[] DPE_;
533	delete[] DPI_;
534	delete[] DPJ_;
535
536	delete[] TBM_;
537	delete[] TBD_;
538	delete[] TBE_;
539	delete[] TBI_;
540	delete[] TBJ_;
541	}
542
543	#if TRACE
544	Log("BestScore=%.6g\n", BestScore);
545	#endif
546	return BestScore;
547	}
548
549	#endif // DOUBLE_AFFINE

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

Context Navigation

source: trunk/GDE/MUSCLE/src/nwdasimple2.cpp

Download in other formats: