+++ /dev/null
-#include "muscle.h"\r
-#include "pwpath.h"\r
-#include "profile.h"\r
-\r
-#if DOUBLE_AFFINE\r
-\r
-#define TRACE 0\r
-\r
-extern bool g_bKeepSimpleDP;\r
-extern SCORE *g_DPM;\r
-extern SCORE *g_DPD;\r
-extern SCORE *g_DPE;\r
-extern SCORE *g_DPI;\r
-extern SCORE *g_DPJ;\r
-extern char *g_TBM;\r
-extern char *g_TBD;\r
-extern char *g_TBE;\r
-extern char *g_TBI;\r
-extern char *g_TBJ;\r
-\r
-static char XlatEdgeType(char c)\r
- {\r
- if ('E' == c)\r
- return 'D';\r
- if ('J' == c)\r
- return 'I';\r
- return c;\r
- }\r
-\r
-static const char *LocalScoreToStr(SCORE s)\r
- {\r
- static char str[16];\r
- if (s < -100000)\r
- return " *";\r
- sprintf(str, "%6.1f", s);\r
- return str;\r
- }\r
-\r
-static void ListDP(const SCORE *DPM_, const ProfPos *PA, const ProfPos *PB,\r
- unsigned uPrefixCountA, unsigned uPrefixCountB)\r
- {\r
- Log(" ");\r
- for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
- {\r
- char c = ' ';\r
- if (uPrefixLengthB > 0)\r
- c = ConsensusChar(PB[uPrefixLengthB - 1]);\r
- Log(" %4u:%c", uPrefixLengthB, c);\r
- }\r
- Log("\n");\r
- for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)\r
- {\r
- char c = ' ';\r
- if (uPrefixLengthA > 0)\r
- c = ConsensusChar(PA[uPrefixLengthA - 1]);\r
- Log("%4u:%c ", uPrefixLengthA, c);\r
- for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
- Log(" %s", LocalScoreToStr(DPM(uPrefixLengthA, uPrefixLengthB)));\r
- Log("\n");\r
- }\r
- }\r
-\r
-static void ListTB(const char *TBM_, const ProfPos *PA, const ProfPos *PB,\r
- unsigned uPrefixCountA, unsigned uPrefixCountB)\r
- {\r
- Log(" ");\r
- for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
- {\r
- char c = ' ';\r
- if (uPrefixLengthB > 0)\r
- c = ConsensusChar(PB[uPrefixLengthB - 1]);\r
- Log(" %4u:%c", uPrefixLengthB, c);\r
- }\r
- Log("\n");\r
- for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)\r
- {\r
- char c = ' ';\r
- if (uPrefixLengthA > 0)\r
- c = ConsensusChar(PA[uPrefixLengthA - 1]);\r
- Log("%4u:%c ", uPrefixLengthA, c);\r
- for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
- Log(" %6c", TBM(uPrefixLengthA, uPrefixLengthB));\r
- Log("\n");\r
- }\r
- }\r
-\r
-static void ListDPM(const SCORE *DPM_, const ProfPos *PA, const ProfPos *PB,\r
- unsigned uPrefixCountA, unsigned uPrefixCountB)\r
- {\r
- Log(" ");\r
- for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
- {\r
- char c = ' ';\r
- if (uPrefixLengthB > 0)\r
- c = ConsensusChar(PB[uPrefixLengthB - 1]);\r
- Log(" %4u:%c", uPrefixLengthB, c);\r
- }\r
- Log("\n");\r
- for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)\r
- {\r
- char c = ' ';\r
- if (uPrefixLengthA > 0)\r
- c = ConsensusChar(PA[uPrefixLengthA - 1]);\r
- Log("%4u:%c ", uPrefixLengthA, c);\r
- for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
- {\r
- SCORE x = (uPrefixLengthA + uPrefixLengthB)*g_scoreGapExtend;\r
- SCORE s = DPM(uPrefixLengthA, uPrefixLengthB) - x;\r
- Log(" %s", LocalScoreToStr(s));\r
- }\r
- Log("\n");\r
- }\r
- }\r
-\r
-extern SCORE ScoreProfPos2(const ProfPos &PP, const ProfPos &PPB);\r
-\r
-SCORE NWDASimple2(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,\r
- unsigned uLengthB, PWPath &Path)\r
- {\r
- assert(uLengthB > 0 && uLengthA > 0);\r
-\r
- const unsigned uPrefixCountA = uLengthA + 1;\r
- const unsigned uPrefixCountB = uLengthB + 1;\r
-\r
-// Allocate DP matrices\r
- const size_t LM = uPrefixCountA*uPrefixCountB;\r
- SCORE *DPM_ = new SCORE[LM];\r
- SCORE *DPD_ = new SCORE[LM];\r
- SCORE *DPE_ = new SCORE[LM];\r
- SCORE *DPI_ = new SCORE[LM];\r
- SCORE *DPJ_ = new SCORE[LM];\r
- SCORE *DPL_ = new SCORE[LM];\r
-\r
- char *TBM_ = new char[LM];\r
- char *TBD_ = new char[LM];\r
- char *TBE_ = new char[LM];\r
- char *TBI_ = new char[LM];\r
- char *TBJ_ = new char[LM];\r
-\r
- memset(DPM_, 0, LM*sizeof(SCORE));\r
- memset(DPD_, 0, LM*sizeof(SCORE));\r
- memset(DPE_, 0, LM*sizeof(SCORE));\r
- memset(DPI_, 0, LM*sizeof(SCORE));\r
- memset(DPJ_, 0, LM*sizeof(SCORE));\r
-\r
-// memset(DPL_, 0, LM*sizeof(SCORE));\r
-\r
- memset(TBM_, '?', LM);\r
- memset(TBD_, '?', LM);\r
- memset(TBE_, '?', LM);\r
- memset(TBI_, '?', LM);\r
- memset(TBJ_, '?', LM);\r
-\r
- DPM(0, 0) = 0;\r
- DPD(0, 0) = MINUS_INFINITY;\r
- DPE(0, 0) = MINUS_INFINITY;\r
- DPI(0, 0) = MINUS_INFINITY;\r
- DPJ(0, 0) = MINUS_INFINITY;\r
-\r
- DPM(1, 0) = MINUS_INFINITY;\r
- DPD(1, 0) = PA[0].m_scoreGapOpen;\r
- DPE(1, 0) = PA[0].m_scoreGapOpen2;\r
- DPI(1, 0) = MINUS_INFINITY;\r
- DPJ(1, 0) = MINUS_INFINITY;\r
-\r
- DPM(0, 1) = MINUS_INFINITY;\r
- DPD(0, 1) = MINUS_INFINITY;\r
- DPE(0, 1) = MINUS_INFINITY;\r
- DPI(0, 1) = PB[0].m_scoreGapOpen;\r
- DPJ(0, 1) = PB[0].m_scoreGapOpen2;\r
-\r
-// Empty prefix of B is special case\r
- for (unsigned uPrefixLengthA = 2; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)\r
- {\r
- // M=LetterA+LetterB, impossible with empty prefix\r
- DPM(uPrefixLengthA, 0) = MINUS_INFINITY;\r
-\r
- // D=LetterA+GapB\r
- DPD(uPrefixLengthA, 0) = DPD(uPrefixLengthA - 1, 0) + g_scoreGapExtend;\r
- TBD(uPrefixLengthA, 0) = 'D';\r
-\r
- DPE(uPrefixLengthA, 0) = DPE(uPrefixLengthA - 1, 0) + g_scoreGapExtend2;\r
- TBE(uPrefixLengthA, 0) = 'E';\r
-\r
- // I=GapA+LetterB, impossible with empty prefix\r
- DPI(uPrefixLengthA, 0) = MINUS_INFINITY;\r
- DPJ(uPrefixLengthA, 0) = MINUS_INFINITY;\r
- }\r
-\r
-// Empty prefix of A is special case\r
- for (unsigned uPrefixLengthB = 2; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
- {\r
- // M=LetterA+LetterB, impossible with empty prefix\r
- DPM(0, uPrefixLengthB) = MINUS_INFINITY;\r
-\r
- // D=LetterA+GapB, impossible with empty prefix\r
- DPD(0, uPrefixLengthB) = MINUS_INFINITY;\r
- DPE(0, uPrefixLengthB) = MINUS_INFINITY;\r
-\r
- // I=GapA+LetterB\r
- DPI(0, uPrefixLengthB) = DPI(0, uPrefixLengthB - 1) + g_scoreGapExtend;\r
- TBI(0, uPrefixLengthB) = 'I';\r
-\r
- DPJ(0, uPrefixLengthB) = DPJ(0, uPrefixLengthB - 1) + g_scoreGapExtend2;\r
- TBJ(0, uPrefixLengthB) = 'J';\r
- }\r
-\r
-// ============\r
-// Main DP loop\r
-// ============\r
- for (unsigned uPrefixLengthB = 1; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
- {\r
- const ProfPos &PPB = PB[uPrefixLengthB - 1];\r
- SCORE scoreGapCloseB;\r
- if (uPrefixLengthB == 1)\r
- scoreGapCloseB = MINUS_INFINITY;\r
- else\r
- scoreGapCloseB = PB[uPrefixLengthB-2].m_scoreGapClose;\r
-\r
- SCORE scoreGapClose2B;\r
- if (uPrefixLengthB == 1)\r
- scoreGapClose2B = MINUS_INFINITY;\r
- else\r
- scoreGapClose2B = PB[uPrefixLengthB-2].m_scoreGapClose2;\r
-\r
- for (unsigned uPrefixLengthA = 1; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)\r
- {\r
- const ProfPos &PPA = PA[uPrefixLengthA - 1];\r
-\r
- {\r
- // Match M=LetterA+LetterB\r
- SCORE scoreLL = ScoreProfPos2(PPA, PPB);\r
- DPL(uPrefixLengthA, uPrefixLengthB) = scoreLL;\r
-\r
- SCORE scoreGapCloseA;\r
- if (uPrefixLengthA == 1)\r
- scoreGapCloseA = MINUS_INFINITY;\r
- else\r
- scoreGapCloseA = PA[uPrefixLengthA-2].m_scoreGapClose;\r
-\r
- SCORE scoreGapClose2A;\r
- if (uPrefixLengthA == 1)\r
- scoreGapClose2A = MINUS_INFINITY;\r
- else\r
- scoreGapClose2A = PA[uPrefixLengthA-2].m_scoreGapClose2;\r
-\r
- SCORE scoreMM = DPM(uPrefixLengthA-1, uPrefixLengthB-1);\r
- SCORE scoreDM = DPD(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseA;\r
- SCORE scoreEM = DPE(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapClose2A;\r
- SCORE scoreIM = DPI(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseB;\r
- SCORE scoreJM = DPJ(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapClose2B;\r
- SCORE scoreBest;\r
- if (scoreMM >= scoreDM && scoreMM >= scoreIM && scoreMM >= scoreEM && scoreMM >= scoreJM)\r
- {\r
- scoreBest = scoreMM;\r
- TBM(uPrefixLengthA, uPrefixLengthB) = 'M';\r
- }\r
- else if (scoreDM >= scoreMM && scoreDM >= scoreIM && scoreDM >= scoreEM && scoreDM >= scoreJM)\r
- {\r
- scoreBest = scoreDM;\r
- TBM(uPrefixLengthA, uPrefixLengthB) = 'D';\r
- }\r
- else if (scoreEM >= scoreMM && scoreEM >= scoreIM && scoreEM >= scoreDM && scoreEM >= scoreJM)\r
- {\r
- scoreBest = scoreEM;\r
- TBM(uPrefixLengthA, uPrefixLengthB) = 'E';\r
- }\r
- else if (scoreIM >= scoreMM && scoreIM >= scoreDM && scoreIM >= scoreEM && scoreIM >= scoreJM)\r
- {\r
- scoreBest = scoreIM;\r
- TBM(uPrefixLengthA, uPrefixLengthB) = 'I';\r
- }\r
- else if (scoreJM >= scoreMM && scoreJM >= scoreDM && scoreJM >= scoreEM && scoreJM >= scoreIM)\r
- {\r
- scoreBest = scoreJM;\r
- TBM(uPrefixLengthA, uPrefixLengthB) = 'J';\r
- }\r
- else\r
- Quit("Max failed (M)");\r
-\r
- DPM(uPrefixLengthA, uPrefixLengthB) = scoreBest + scoreLL;\r
- }\r
-\r
- {\r
- // Delete D=LetterA+GapB\r
- SCORE scoreMD = DPM(uPrefixLengthA-1, uPrefixLengthB) +\r
- PA[uPrefixLengthA-1].m_scoreGapOpen;\r
- SCORE scoreDD = DPD(uPrefixLengthA-1, uPrefixLengthB) +\r
- g_scoreGapExtend;\r
-\r
- SCORE scoreBest;\r
- if (scoreMD >= scoreDD)\r
- {\r
- scoreBest = scoreMD;\r
- TBD(uPrefixLengthA, uPrefixLengthB) = 'M';\r
- }\r
- else\r
- {\r
- assert(scoreDD >= scoreMD);\r
- scoreBest = scoreDD;\r
- TBD(uPrefixLengthA, uPrefixLengthB) = 'D';\r
- }\r
- DPD(uPrefixLengthA, uPrefixLengthB) = scoreBest;\r
- }\r
-\r
- {\r
- // Delete E=LetterA+GapB\r
- SCORE scoreME = DPM(uPrefixLengthA-1, uPrefixLengthB) +\r
- PA[uPrefixLengthA-1].m_scoreGapOpen2;\r
- SCORE scoreEE = DPE(uPrefixLengthA-1, uPrefixLengthB) +\r
- g_scoreGapExtend2;\r
-\r
- SCORE scoreBest;\r
- if (scoreME >= scoreEE)\r
- {\r
- scoreBest = scoreME;\r
- TBE(uPrefixLengthA, uPrefixLengthB) = 'M';\r
- }\r
- else\r
- {\r
- assert(scoreEE >= scoreME);\r
- scoreBest = scoreEE;\r
- TBE(uPrefixLengthA, uPrefixLengthB) = 'E';\r
- }\r
- DPE(uPrefixLengthA, uPrefixLengthB) = scoreBest;\r
- }\r
-\r
- // Insert I=GapA+LetterB\r
- {\r
- SCORE scoreMI = DPM(uPrefixLengthA, uPrefixLengthB-1) +\r
- PB[uPrefixLengthB-1].m_scoreGapOpen;\r
- SCORE scoreII = DPI(uPrefixLengthA, uPrefixLengthB-1) +\r
- g_scoreGapExtend;\r
-\r
- SCORE scoreBest;\r
- if (scoreMI >= scoreII)\r
- {\r
- scoreBest = scoreMI;\r
- TBI(uPrefixLengthA, uPrefixLengthB) = 'M';\r
- }\r
- else \r
- {\r
- assert(scoreII > scoreMI);\r
- scoreBest = scoreII;\r
- TBI(uPrefixLengthA, uPrefixLengthB) = 'I';\r
- }\r
- DPI(uPrefixLengthA, uPrefixLengthB) = scoreBest;\r
- }\r
-\r
- // Insert J=GapA+LetterB\r
- {\r
- SCORE scoreMJ = DPM(uPrefixLengthA, uPrefixLengthB-1) +\r
- PB[uPrefixLengthB-1].m_scoreGapOpen2;\r
- SCORE scoreJJ = DPJ(uPrefixLengthA, uPrefixLengthB-1) +\r
- g_scoreGapExtend2;\r
-\r
- SCORE scoreBest;\r
- if (scoreMJ > scoreJJ)\r
- {\r
- scoreBest = scoreMJ;\r
- TBJ(uPrefixLengthA, uPrefixLengthB) = 'M';\r
- }\r
- else \r
- {\r
- assert(scoreJJ >= scoreMJ);\r
- scoreBest = scoreJJ;\r
- TBJ(uPrefixLengthA, uPrefixLengthB) = 'J';\r
- }\r
- DPJ(uPrefixLengthA, uPrefixLengthB) = scoreBest;\r
- }\r
- }\r
- }\r
-\r
-// Special case: close gaps at end of alignment\r
- DPD(uLengthA, uLengthB) += PA[uLengthA-1].m_scoreGapClose;\r
- DPE(uLengthA, uLengthB) += PA[uLengthA-1].m_scoreGapClose2;\r
-\r
- DPI(uLengthA, uLengthB) += PB[uLengthB-1].m_scoreGapClose;\r
- DPJ(uLengthA, uLengthB) += PB[uLengthB-1].m_scoreGapClose2;\r
-\r
-#if TRACE\r
- Log("DPL:\n");\r
- ListDP(DPL_, PA, PB, uPrefixCountA, uPrefixCountB);\r
- Log("DPM:\n");\r
- ListDP(DPM_, PA, PB, uPrefixCountA, uPrefixCountB);\r
- Log("DPD:\n");\r
- ListDP(DPD_, PA, PB, uPrefixCountA, uPrefixCountB);\r
- Log("DPE:\n");\r
- ListDP(DPE_, PA, PB, uPrefixCountA, uPrefixCountB);\r
- Log("DPI:\n");\r
- ListDP(DPI_, PA, PB, uPrefixCountA, uPrefixCountB);\r
- Log("DPJ:\n");\r
- ListDP(DPJ_, PA, PB, uPrefixCountA, uPrefixCountB);\r
- Log("TBM:\n");\r
- ListTB(TBM_, PA, PB, uPrefixCountA, uPrefixCountB);\r
- Log("TBD:\n");\r
- ListTB(TBD_, PA, PB, uPrefixCountA, uPrefixCountB);\r
- Log("TBE:\n");\r
- ListTB(TBE_, PA, PB, uPrefixCountA, uPrefixCountB);\r
- Log("TBI:\n");\r
- ListTB(TBI_, PA, PB, uPrefixCountA, uPrefixCountB);\r
- Log("TBJ:\n");\r
- ListTB(TBJ_, PA, PB, uPrefixCountA, uPrefixCountB);\r
-#endif\r
-\r
-// ==========\r
-// Trace-back\r
-// ==========\r
-\r
- Path.Clear();\r
-\r
-// Find last edge\r
- char cEdgeType = '?';\r
- SCORE BestScore = MINUS_INFINITY;\r
- SCORE M = DPM(uLengthA, uLengthB);\r
- SCORE D = DPD(uLengthA, uLengthB);\r
- SCORE E = DPE(uLengthA, uLengthB);\r
- SCORE I = DPI(uLengthA, uLengthB);\r
- SCORE J = DPJ(uLengthA, uLengthB);\r
-\r
- if (M >= D && M >= E && M >= I && M >= J)\r
- {\r
- cEdgeType = 'M';\r
- BestScore = M;\r
- }\r
- else if (D >= M && D >= E && D >= I && D >= J)\r
- {\r
- cEdgeType = 'D';\r
- BestScore = D;\r
- }\r
- else if (E >= M && E >= D && E >= I && E >= J)\r
- {\r
- cEdgeType = 'E';\r
- BestScore = E;\r
- }\r
- else if (I >= M && I >= D && I >= E && I >= J)\r
- {\r
- cEdgeType = 'I';\r
- BestScore = I;\r
- }\r
- else if (J >= M && J >= D && J >= E && J >= I)\r
- {\r
- cEdgeType = 'J';\r
- BestScore = J;\r
- }\r
- else\r
- Quit("Bad max");\r
-\r
- unsigned PLA = uLengthA;\r
- unsigned PLB = uLengthB;\r
- unsigned ECount = 0;\r
- unsigned JCount = 0;\r
- for (;;)\r
- {\r
-#if TRACE\r
- Log("TraceBack: %c%u.%u\n", cEdgeType, PLA, PLB);\r
-#endif\r
- PWEdge Edge;\r
- Edge.cType = XlatEdgeType(cEdgeType);\r
- Edge.uPrefixLengthA = PLA;\r
- Edge.uPrefixLengthB = PLB;\r
- Path.PrependEdge(Edge);\r
-\r
- switch (cEdgeType)\r
- {\r
- case 'M':\r
- assert(PLA > 0);\r
- assert(PLB > 0);\r
- cEdgeType = TBM(PLA, PLB);\r
- --PLA;\r
- --PLB;\r
- break;\r
-\r
- case 'D':\r
- assert(PLA > 0);\r
- cEdgeType = TBD(PLA, PLB);\r
- --PLA;\r
- break;\r
-\r
- case 'E':\r
- ++ECount;\r
- assert(PLA > 0);\r
- cEdgeType = TBE(PLA, PLB);\r
- --PLA;\r
- break;\r
-\r
- case 'I':\r
- assert(PLB > 0);\r
- cEdgeType = TBI(PLA, PLB);\r
- --PLB;\r
- break;\r
- \r
- case 'J':\r
- ++JCount;\r
- assert(PLB > 0);\r
- cEdgeType = TBJ(PLA, PLB);\r
- --PLB;\r
- break;\r
-\r
- default:\r
- Quit("Invalid edge %c", cEdgeType);\r
- }\r
- if (0 == PLA && 0 == PLB)\r
- break;\r
- }\r
- //if (ECount > 0 || JCount > 0)\r
- // fprintf(stderr, "E=%d J=%d\n", ECount, JCount);\r
- Path.Validate();\r
- if (Path.GetMatchCount() + Path.GetDeleteCount() != uLengthA)\r
- Quit("Path count A");\r
- if (Path.GetMatchCount() + Path.GetInsertCount() != uLengthB)\r
- Quit("Path count B");\r
-\r
- if (g_bKeepSimpleDP)\r
- {\r
- g_DPM = DPM_;\r
- g_DPD = DPD_;\r
- g_DPE = DPE_;\r
- g_DPI = DPI_;\r
- g_DPJ = DPJ_;\r
-\r
- g_TBM = TBM_;\r
- g_TBD = TBD_;\r
- g_TBE = TBE_;\r
- g_TBI = TBI_;\r
- g_TBJ = TBJ_;\r
- }\r
- else\r
- {\r
- delete[] DPM_;\r
- delete[] DPD_;\r
- delete[] DPE_;\r
- delete[] DPI_;\r
- delete[] DPJ_;\r
-\r
- delete[] TBM_;\r
- delete[] TBD_;\r
- delete[] TBE_;\r
- delete[] TBI_;\r
- delete[] TBJ_;\r
- }\r
-\r
-#if TRACE\r
- Log("BestScore=%.6g\n", BestScore);\r
-#endif\r
- return BestScore;\r
- }\r
-\r
-#endif // DOUBLE_AFFINE\r