+++ /dev/null
-#include "muscle.h"\r
-#include <math.h>\r
-#include "pwpath.h"\r
-#include "profile.h"\r
-#include <stdio.h>\r
-\r
-#define TRACE 0\r
-\r
-bool g_bKeepSimpleDP;\r
-SCORE *g_DPM;\r
-SCORE *g_DPD;\r
-SCORE *g_DPE;\r
-SCORE *g_DPI;\r
-SCORE *g_DPJ;\r
-char *g_TBM;\r
-char *g_TBD;\r
-char *g_TBE;\r
-char *g_TBI;\r
-char *g_TBJ;\r
-\r
-#if DOUBLE_AFFINE\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 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 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
-SCORE NWDASimple(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 *DPL_ = new SCORE[LM];\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
-\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(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
- TBD(1, 0) = 'D';\r
- TBE(1, 0) = 'E';\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
- TBI(0, 1) = 'I';\r
- TBJ(0, 1) = 'J';\r
-\r
-// Empty prefix of B is special case\r
- for (unsigned uPrefixLengthA = 2; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)\r
- {\r
- DPM(uPrefixLengthA, 0) = MINUS_INFINITY;\r
-\r
- DPD(uPrefixLengthA, 0) = DPD(uPrefixLengthA - 1, 0) + g_scoreGapExtend;\r
- DPE(uPrefixLengthA, 0) = DPE(uPrefixLengthA - 1, 0) + g_scoreGapExtend2;\r
-\r
- TBD(uPrefixLengthA, 0) = 'D';\r
- TBE(uPrefixLengthA, 0) = 'E';\r
-\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
- DPM(0, uPrefixLengthB) = MINUS_INFINITY;\r
-\r
- DPD(0, uPrefixLengthB) = MINUS_INFINITY;\r
- DPE(0, uPrefixLengthB) = MINUS_INFINITY;\r
-\r
- DPI(0, uPrefixLengthB) = DPI(0, uPrefixLengthB - 1) + g_scoreGapExtend;\r
- DPJ(0, uPrefixLengthB) = DPJ(0, uPrefixLengthB - 1) + g_scoreGapExtend2;\r
-\r
- TBI(0, uPrefixLengthB) = 'I';\r
- TBJ(0, uPrefixLengthB) = 'J';\r
- }\r
-\r
-// Special case to agree with NWFast, no D-I transitions so...\r
- DPD(uLengthA, 0) = MINUS_INFINITY;\r
- DPE(uLengthA, 0) = MINUS_INFINITY;\r
-// DPI(0, uLengthB) = MINUS_INFINITY;\r
-// DPJ(0, uLengthB) = MINUS_INFINITY;\r
-\r
-// ============\r
-// Main DP loop\r
-// ============\r
- SCORE scoreGapCloseB = MINUS_INFINITY;\r
- SCORE scoreGapClose2B = MINUS_INFINITY;\r
- for (unsigned uPrefixLengthB = 1; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
- {\r
- const ProfPos &PPB = PB[uPrefixLengthB - 1];\r
-\r
- SCORE scoreGapCloseA = MINUS_INFINITY;\r
- SCORE scoreGapClose2A = MINUS_INFINITY;\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 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
-\r
- SCORE scoreBest;\r
- if (scoreMM >= scoreDM && scoreMM >= scoreEM && scoreMM >= scoreIM && scoreMM >= scoreJM)\r
- {\r
- scoreBest = scoreMM;\r
- TBM(uPrefixLengthA, uPrefixLengthB) = 'M';\r
- }\r
- else if (scoreDM >= scoreMM && scoreDM >= scoreEM && scoreDM >= scoreIM && scoreDM >= scoreJM)\r
- {\r
- scoreBest = scoreDM;\r
- TBM(uPrefixLengthA, uPrefixLengthB) = 'D';\r
- }\r
- else if (scoreEM >= scoreMM && scoreEM >= scoreDM && scoreEM >= scoreIM && 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\r
- {\r
- assert(scoreJM >= scoreMM && scoreJM >= scoreDM && scoreJM >= scoreEM && scoreJM >= scoreIM);\r
- scoreBest = scoreJM;\r
- TBM(uPrefixLengthA, uPrefixLengthB) = 'J';\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) + 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) + 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) + 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) + 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
- scoreGapCloseA = PPA.m_scoreGapClose;\r
- scoreGapClose2A = PPA.m_scoreGapClose2;\r
- }\r
- scoreGapCloseB = PPB.m_scoreGapClose;\r
- scoreGapClose2B = PPB.m_scoreGapClose2;\r
- }\r
-\r
-#if TRACE\r
- Log("\n");\r
- Log("DA Simple DPL:\n");\r
- ListDP(DPL_, PA, PB, uPrefixCountA, uPrefixCountB);\r
- Log("\n");\r
- Log("DA Simple DPM:\n");\r
- ListDP(DPM_, PA, PB, uPrefixCountA, uPrefixCountB);\r
- Log("\n");\r
- Log("DA Simple DPD:\n");\r
- ListDP(DPD_, PA, PB, uPrefixCountA, uPrefixCountB);\r
- Log("\n");\r
- Log("DA Simple DPE:\n");\r
- ListDP(DPE_, PA, PB, uPrefixCountA, uPrefixCountB);\r
- Log("\n");\r
- Log("DA Simple DPI:\n");\r
- ListDP(DPI_, PA, PB, uPrefixCountA, uPrefixCountB);\r
- Log("\n");\r
- Log("DA Simple DPJ:\n");\r
- ListDP(DPJ_, PA, PB, uPrefixCountA, uPrefixCountB);\r
- Log("\n");\r
- Log("DA Simple TBM:\n");\r
- ListTB(TBM_, PA, PB, uPrefixCountA, uPrefixCountB);\r
- Log("\n");\r
- Log("DA Simple TBD:\n");\r
- ListTB(TBD_, PA, PB, uPrefixCountA, uPrefixCountB);\r
- Log("\n");\r
- Log("DA Simple TBE:\n");\r
- ListTB(TBE_, PA, PB, uPrefixCountA, uPrefixCountB);\r
- Log("\n");\r
- Log("DA Simple TBI:\n");\r
- ListTB(TBI_, PA, PB, uPrefixCountA, uPrefixCountB);\r
- Log("\n");\r
- Log("DA Simple TBJ:\n");\r
- ListTB(TBJ_, PA, PB, uPrefixCountA, uPrefixCountB);\r
-#endif\r
-\r
-// Trace-back\r
-// ==========\r
- Path.Clear();\r
-\r
-// Find last edge\r
- SCORE M = DPM(uLengthA, uLengthB);\r
- SCORE D = DPD(uLengthA, uLengthB) + PA[uLengthA-1].m_scoreGapClose;\r
- SCORE E = DPE(uLengthA, uLengthB) + PA[uLengthA-1].m_scoreGapClose2;\r
- SCORE I = DPI(uLengthA, uLengthB) + PB[uLengthB-1].m_scoreGapClose;\r
- SCORE J = DPJ(uLengthA, uLengthB) + PB[uLengthB-1].m_scoreGapClose2;\r
- char cEdgeType = '?';\r
-\r
- SCORE BestScore = M;\r
- cEdgeType = 'M';\r
- if (D > BestScore)\r
- {\r
- cEdgeType = 'D';\r
- BestScore = D;\r
- }\r
- if (E > BestScore)\r
- {\r
- cEdgeType = 'E';\r
- BestScore = E;\r
- }\r
- if (I > BestScore)\r
- {\r
- cEdgeType = 'I';\r
- BestScore = I;\r
- }\r
- if (J > BestScore)\r
- {\r
- cEdgeType = 'J';\r
- BestScore = J;\r
- }\r
-\r
-#if TRACE\r
- Log("DA Simple: MAB=%.4g DAB=%.4g EAB=%.4g IAB=%.4g JAB=%.4g best=%c\n",\r
- M, D, E, I, J, cEdgeType);\r
-#endif\r
-\r
- unsigned PLA = uLengthA;\r
- unsigned PLB = uLengthB;\r
- for (;;)\r
- {\r
- PWEdge Edge;\r
- Edge.cType = XlatEdgeType(cEdgeType);\r
- Edge.uPrefixLengthA = PLA;\r
- Edge.uPrefixLengthB = PLB;\r
-#if TRACE\r
- Log("Prepend %c%d.%d\n", Edge.cType, PLA, PLB);\r
-#endif\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
- 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
- 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
- Path.Validate();\r
-\r
-// SCORE Score = TraceBack(PA, uLengthA, PB, uLengthB, DPM_, DPD_, DPI_, Path);\r
-\r
-#if TRACE\r
- SCORE scorePath = FastScorePath2(PA, uLengthA, PB, uLengthB, Path);\r
- Path.LogMe();\r
- Log("Score = %s Path = %s\n", LocalScoreToStr(BestScore), LocalScoreToStr(scorePath));\r
-#endif\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
- return BestScore;\r
- }\r
-\r
-#endif // DOUBLE_AFFINE\r