+++ /dev/null
-#include "muscle.h"\r
-#include "profile.h"\r
-#include "pwpath.h"\r
-#include <math.h>\r
-\r
-#define TRACE 0\r
-\r
-#define EQ(a, b) (fabs(a-b) < 0.1)\r
-\r
-SCORE TraceBack(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,\r
- unsigned uLengthB, const SCORE *DPM_, const SCORE *DPD_, const SCORE *DPI_,\r
- PWPath &Path)\r
- {\r
-#if TRACE\r
- Log("\n");\r
- Log("TraceBack LengthA=%u LengthB=%u\n", uLengthA, uLengthB);\r
-#endif\r
- assert(uLengthB > 0 && uLengthA > 0);\r
-\r
- const unsigned uPrefixCountA = uLengthA + 1;\r
- const unsigned uPrefixCountB = uLengthB + 1;\r
-\r
- Path.Clear();\r
-\r
- unsigned uPrefixLengthA = uLengthA;\r
- unsigned uPrefixLengthB = uLengthB;\r
-\r
- const SCORE scoreM = DPM(uPrefixLengthA, uPrefixLengthB);\r
- SCORE scoreD = DPD(uPrefixLengthA, uPrefixLengthB);\r
- SCORE scoreI = DPI(uPrefixLengthA, uPrefixLengthB);\r
-\r
- const ProfPos &LastPPA = PA[uLengthA - 1];\r
- const ProfPos &LastPPB = PB[uLengthB - 1];\r
-\r
- scoreD += LastPPA.m_scoreGapClose;\r
- scoreI += LastPPB.m_scoreGapClose;\r
-\r
- char cEdgeType = cInsane;\r
- SCORE scoreMax;\r
- if (scoreM >= scoreD && scoreM >= scoreI)\r
- {\r
- scoreMax = scoreM;\r
- cEdgeType = 'M';\r
- }\r
- else if (scoreD >= scoreM && scoreD >= scoreI)\r
- {\r
- scoreMax = scoreD;\r
- cEdgeType = 'D';\r
- }\r
- else\r
- {\r
- assert(scoreI >= scoreM && scoreI >= scoreD);\r
- scoreMax = scoreI;\r
- cEdgeType = 'I';\r
- }\r
-\r
- for (;;)\r
- {\r
- if ('S' == cEdgeType)\r
- break;\r
-\r
- PWEdge Edge;\r
- Edge.cType = cEdgeType;\r
- Edge.uPrefixLengthA = uPrefixLengthA;\r
- Edge.uPrefixLengthB = uPrefixLengthB;\r
- Path.PrependEdge(Edge);\r
-\r
- char cPrevEdgeType;\r
- unsigned uPrevPrefixLengthA = uPrefixLengthA;\r
- unsigned uPrevPrefixLengthB = uPrefixLengthB;\r
-\r
- switch (cEdgeType)\r
- {\r
- case 'M':\r
- {\r
- assert(uPrefixLengthA > 0);\r
- assert(uPrefixLengthB > 0);\r
- const ProfPos &PPA = PA[uPrefixLengthA - 1];\r
- const ProfPos &PPB = PB[uPrefixLengthB - 1];\r
-\r
- const SCORE Score = DPM(uPrefixLengthA, uPrefixLengthB);\r
- const SCORE scoreMatch = ScoreProfPos2(PPA, PPB);\r
-\r
- SCORE scoreSM;\r
- if (1 == uPrefixLengthA && 1 == uPrefixLengthB)\r
- scoreSM = scoreMatch;\r
- else\r
- scoreSM = MINUS_INFINITY;\r
-\r
- SCORE scoreMM = MINUS_INFINITY;\r
- SCORE scoreDM = MINUS_INFINITY;\r
- SCORE scoreIM = MINUS_INFINITY;\r
- if (uPrefixLengthA > 1 && uPrefixLengthB > 1)\r
- scoreMM = DPM(uPrefixLengthA-1, uPrefixLengthB-1) + scoreMatch;\r
- if (uPrefixLengthA > 1)\r
- {\r
- SCORE scoreTransDM = PA[uPrefixLengthA-2].m_scoreGapClose;\r
- scoreDM = DPD(uPrefixLengthA-1, uPrefixLengthB-1) + scoreTransDM + scoreMatch;\r
- }\r
- if (uPrefixLengthB > 1)\r
- {\r
- SCORE scoreTransIM = PB[uPrefixLengthB-2].m_scoreGapClose;\r
- scoreIM = DPI(uPrefixLengthA-1, uPrefixLengthB-1) + scoreTransIM + scoreMatch;\r
- }\r
-\r
- if (EQ(scoreMM, Score))\r
- cPrevEdgeType = 'M';\r
- else if (EQ(scoreDM, Score))\r
- cPrevEdgeType = 'D';\r
- else if (EQ(scoreIM, Score))\r
- cPrevEdgeType = 'I';\r
- else if (EQ(scoreSM, Score))\r
- cPrevEdgeType = 'S';\r
- else\r
- Quit("TraceBack: failed to match M score=%g M=%g D=%g I=%g S=%g",\r
- Score, scoreMM, scoreDM, scoreIM, scoreSM);\r
-\r
- --uPrevPrefixLengthA;\r
- --uPrevPrefixLengthB;\r
- break;\r
- }\r
-\r
- case 'D':\r
- {\r
- assert(uPrefixLengthA > 0);\r
- const SCORE Score = DPD(uPrefixLengthA, uPrefixLengthB);\r
-\r
- SCORE scoreMD = MINUS_INFINITY;\r
- SCORE scoreDD = MINUS_INFINITY;\r
- SCORE scoreSD = MINUS_INFINITY;\r
- if (uPrefixLengthB == 0)\r
- {\r
- if (uPrefixLengthA == 1)\r
- scoreSD = PA[0].m_scoreGapOpen;\r
- else\r
- scoreSD = DPD(uPrefixLengthA - 1, 0);\r
- }\r
- if (uPrefixLengthA > 1)\r
- {\r
- const ProfPos &PPA = PA[uPrefixLengthA - 1];\r
- SCORE scoreTransMD = PPA.m_scoreGapOpen;\r
- scoreMD = DPM(uPrefixLengthA-1, uPrefixLengthB) + scoreTransMD;\r
- scoreDD = DPD(uPrefixLengthA-1, uPrefixLengthB);\r
- }\r
-\r
- if (EQ(Score, scoreMD))\r
- cPrevEdgeType = 'M';\r
- else if (EQ(Score, scoreDD))\r
- cPrevEdgeType = 'D';\r
- else if (EQ(Score, scoreSD))\r
- cPrevEdgeType = 'S';\r
- else\r
- Quit("TraceBack: failed to match D");\r
-\r
- --uPrevPrefixLengthA;\r
- break;\r
- }\r
-\r
- case 'I':\r
- {\r
- assert(uPrefixLengthB > 0);\r
- const SCORE Score = DPI(uPrefixLengthA, uPrefixLengthB);\r
-\r
- SCORE scoreMI = MINUS_INFINITY;\r
- SCORE scoreII = MINUS_INFINITY;\r
- SCORE scoreSI = MINUS_INFINITY;\r
- if (uPrefixLengthA == 0)\r
- {\r
- if (uPrefixLengthB == 1)\r
- scoreSI = PB[0].m_scoreGapOpen;\r
- else\r
- scoreSI = DPI(0, uPrefixLengthB - 1);\r
- }\r
- if (uPrefixLengthB > 1)\r
- {\r
- const ProfPos &PPB = PB[uPrefixLengthB - 1];\r
- SCORE scoreTransMI = PPB.m_scoreGapOpen;\r
- scoreMI = DPM(uPrefixLengthA, uPrefixLengthB-1) + scoreTransMI;\r
- scoreII = DPI(uPrefixLengthA, uPrefixLengthB-1);\r
- }\r
-\r
- if (EQ(Score, scoreMI))\r
- cPrevEdgeType = 'M';\r
- else if (EQ(Score, scoreII))\r
- cPrevEdgeType = 'I';\r
- else if (EQ(Score, scoreSI))\r
- cPrevEdgeType = 'S';\r
- else\r
- Quit("TraceBack: failed to match I");\r
-\r
- --uPrevPrefixLengthB;\r
- break;\r
- }\r
-\r
- default:\r
- assert(false);\r
- }\r
-#if TRACE\r
- Log("Edge %c%c%u.%u", cPrevEdgeType, cEdgeType, uPrefixLengthA, uPrefixLengthB);\r
- Log("\n");\r
-#endif\r
- cEdgeType = cPrevEdgeType;\r
- uPrefixLengthA = uPrevPrefixLengthA;\r
- uPrefixLengthB = uPrevPrefixLengthB;\r
- }\r
-\r
- return scoreMax;\r
- }\r
git://source.jalview.org / jabaws.git / blobdiff