source: branches/lib/GDE/MUSCLE/src/nwdasimple.cpp

Last change on this file was 10390, checked in by aboeckma, 12 years ago

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