--- /dev/null
+#include "muscle.h"\r
+#include "msa.h"\r
+#include "pwpath.h"\r
+#include "profile.h"\r
+\r
+#define TRACE 0\r
+\r
+static void LogPP(const ProfPos &PP)\r
+ {\r
+ Log("ResidueGroup %u\n", PP.m_uResidueGroup);\r
+ Log("AllGaps %d\n", PP.m_bAllGaps);\r
+ Log("Occ %.3g\n", PP.m_fOcc);\r
+ Log("LL=%.3g LG=%.3g GL=%.3g GG=%.3g\n", PP.m_LL, PP.m_LG, PP.m_GL, PP.m_GG);\r
+ Log("Freqs ");\r
+ for (unsigned i = 0; i < 20; ++i)\r
+ if (PP.m_fcCounts[i] > 0)\r
+ Log("%c=%.3g ", LetterToChar(i), PP.m_fcCounts[i]);\r
+ Log("\n");\r
+ }\r
+\r
+static void AssertProfPosEq(const ProfPos *PA, const ProfPos *PB, unsigned i)\r
+ {\r
+ const ProfPos &PPA = PA[i];\r
+ const ProfPos &PPB = PB[i];\r
+#define eq(x) if (PPA.m_##x != PPB.m_##x) { LogPP(PPA); LogPP(PPB); Quit("AssertProfPosEq." #x); }\r
+#define be(x) if (!BTEq(PPA.m_##x, PPB.m_##x)) { LogPP(PPA); LogPP(PPB); Quit("AssertProfPosEq." #x); }\r
+ eq(bAllGaps)\r
+ eq(uResidueGroup)\r
+\r
+ be(LL)\r
+ be(LG)\r
+ be(GL)\r
+ be(GG)\r
+ be(fOcc)\r
+ be(scoreGapOpen)\r
+ be(scoreGapClose)\r
+\r
+ for (unsigned j = 0; j < 20; ++j)\r
+ {\r
+#define eqj(x) if (PPA.m_##x != PPB.m_##x) Quit("AssertProfPosEq j=%u " #x, j);\r
+#define bej(x) if (!BTEq(PPA.m_##x, PPB.m_##x)) Quit("AssertProfPosEq j=%u " #x, j);\r
+ bej(fcCounts[j]);\r
+// eqj(uSortOrder[j]) // may differ due to ties, don't check?\r
+ bej(AAScores[j])\r
+#undef eqj\r
+#undef bej\r
+ }\r
+#undef eq\r
+#undef be\r
+ }\r
+\r
+void AssertProfsEq(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,\r
+ unsigned uLengthB)\r
+ {\r
+ if (uLengthA != uLengthB)\r
+ Quit("AssertProfsEq: lengths differ %u %u", uLengthA, uLengthB);\r
+ for (unsigned i = 0; i < uLengthB; ++i)\r
+ AssertProfPosEq(PA, PB, i);\r
+ }\r
+\r
+#if DEBUG\r
+static void ValidateProf(const ProfPos *Prof, unsigned uLength)\r
+ {\r
+ for (unsigned i = 0; i < uLength; ++i)\r
+ {\r
+ const ProfPos &PP = Prof[i];\r
+\r
+ FCOUNT s1 = PP.m_LL + PP.m_LG + PP.m_GL + PP.m_GG;\r
+ assert(BTEq(s1, 1.0));\r
+\r
+ if (i > 0)\r
+ {\r
+ const ProfPos &PPPrev = Prof[i-1];\r
+ FCOUNT s2 = PPPrev.m_LL + PPPrev.m_GL;\r
+ FCOUNT s3 = PP.m_LL + PP.m_LG;\r
+ assert(BTEq(s2, s3));\r
+ }\r
+ if (i < uLength - 1)\r
+ {\r
+ const ProfPos &PPNext = Prof[i+1];\r
+ FCOUNT s4 = PP.m_LL + PP.m_GL;\r
+ FCOUNT s5 = PPNext.m_LL + PPNext.m_LG;\r
+ assert(BTEq(s4, s5));\r
+ }\r
+ }\r
+ }\r
+#else\r
+#define ValidateProf(Prof, Length) /* empty */\r
+#endif\r
+\r
+static void ScoresFromFreqsPos(ProfPos *Prof, unsigned uLength, unsigned uPos)\r
+ {\r
+ ProfPos &PP = Prof[uPos];\r
+ SortCounts(PP.m_fcCounts, PP.m_uSortOrder);\r
+ PP.m_uResidueGroup = ResidueGroupFromFCounts(PP.m_fcCounts);\r
+\r
+// "Occupancy"\r
+ PP.m_fOcc = PP.m_LL + PP.m_GL;\r
+\r
+// Frequency of gap-opens in this position (i)\r
+// Gap open = letter in i-1 and gap in i\r
+// = iff LG in i\r
+ FCOUNT fcOpen = PP.m_LG;\r
+\r
+// Frequency of gap-closes in this position\r
+// Gap close = gap in i and letter in i+1\r
+// = iff GL in i+1\r
+ FCOUNT fcClose;\r
+ if (uPos + 1 < uLength)\r
+ fcClose = Prof[uPos + 1].m_GL;\r
+ else\r
+ fcClose = PP.m_GG + PP.m_LG;\r
+\r
+ PP.m_scoreGapOpen = (SCORE) ((1.0 - fcOpen)*g_scoreGapOpen/2.0);\r
+ PP.m_scoreGapClose = (SCORE) ((1.0 - fcClose)*g_scoreGapOpen/2.0);\r
+#if DOUBLE_AFFINE\r
+ PP.m_scoreGapOpen2 = (SCORE) ((1.0 - fcOpen)*g_scoreGapOpen2/2.0);\r
+ PP.m_scoreGapClose2 = (SCORE) ((1.0 - fcClose)*g_scoreGapOpen2/2.0);\r
+#endif\r
+\r
+ for (unsigned i = 0; i < g_AlphaSize; ++i)\r
+ {\r
+ SCORE scoreSum = 0;\r
+ for (unsigned j = 0; j < g_AlphaSize; ++j)\r
+ scoreSum += PP.m_fcCounts[j]*(*g_ptrScoreMatrix)[i][j];\r
+ PP.m_AAScores[i] = scoreSum;\r
+ }\r
+ }\r
+\r
+void ProfScoresFromFreqs(ProfPos *Prof, unsigned uLength)\r
+ {\r
+ for (unsigned i = 0; i < uLength; ++i)\r
+ ScoresFromFreqsPos(Prof, uLength, i);\r
+ }\r
+\r
+static void AppendDelete(const MSA &msaA, unsigned &uColIndexA,\r
+ unsigned uSeqCountA, unsigned uSeqCountB, MSA &msaCombined,\r
+ unsigned &uColIndexCombined)\r
+ {\r
+#if TRACE\r
+ Log("AppendDelete ColIxA=%u ColIxCmb=%u\n",\r
+ uColIndexA, uColIndexCombined);\r
+#endif\r
+ for (unsigned uSeqIndexA = 0; uSeqIndexA < uSeqCountA; ++uSeqIndexA)\r
+ {\r
+ char c = msaA.GetChar(uSeqIndexA, uColIndexA);\r
+ msaCombined.SetChar(uSeqIndexA, uColIndexCombined, c);\r
+ }\r
+\r
+ for (unsigned uSeqIndexB = 0; uSeqIndexB < uSeqCountB; ++uSeqIndexB)\r
+ msaCombined.SetChar(uSeqCountA + uSeqIndexB, uColIndexCombined, '-');\r
+\r
+ ++uColIndexCombined;\r
+ ++uColIndexA;\r
+ }\r
+\r
+static void AppendInsert(const MSA &msaB, unsigned &uColIndexB,\r
+ unsigned uSeqCountA, unsigned uSeqCountB, MSA &msaCombined,\r
+ unsigned &uColIndexCombined)\r
+ {\r
+#if TRACE\r
+ Log("AppendInsert ColIxB=%u ColIxCmb=%u\n",\r
+ uColIndexB, uColIndexCombined);\r
+#endif\r
+ for (unsigned uSeqIndexA = 0; uSeqIndexA < uSeqCountA; ++uSeqIndexA)\r
+ msaCombined.SetChar(uSeqIndexA, uColIndexCombined, '-');\r
+\r
+ for (unsigned uSeqIndexB = 0; uSeqIndexB < uSeqCountB; ++uSeqIndexB)\r
+ {\r
+ char c = msaB.GetChar(uSeqIndexB, uColIndexB);\r
+ msaCombined.SetChar(uSeqCountA + uSeqIndexB, uColIndexCombined, c);\r
+ }\r
+\r
+ ++uColIndexCombined;\r
+ ++uColIndexB;\r
+ }\r
+\r
+static void AppendTplInserts(const MSA &msaA, unsigned &uColIndexA, unsigned uColCountA,\r
+ const MSA &msaB, unsigned &uColIndexB, unsigned uColCountB, unsigned uSeqCountA,\r
+ unsigned uSeqCountB, MSA &msaCombined, unsigned &uColIndexCombined)\r
+ {\r
+#if TRACE\r
+ Log("AppendTplInserts ColIxA=%u ColIxB=%u ColIxCmb=%u\n",\r
+ uColIndexA, uColIndexB, uColIndexCombined);\r
+#endif\r
+ const unsigned uLengthA = msaA.GetColCount();\r
+ const unsigned uLengthB = msaB.GetColCount();\r
+\r
+ unsigned uNewColCount = uColCountA;\r
+ if (uColCountB > uNewColCount)\r
+ uNewColCount = uColCountB;\r
+\r
+ for (unsigned n = 0; n < uColCountA; ++n)\r
+ {\r
+ for (unsigned uSeqIndexA = 0; uSeqIndexA < uSeqCountA; ++uSeqIndexA)\r
+ {\r
+ char c = msaA.GetChar(uSeqIndexA, uColIndexA + n);\r
+ c = UnalignChar(c);\r
+ msaCombined.SetChar(uSeqIndexA, uColIndexCombined + n, c);\r
+ }\r
+ }\r
+ for (unsigned n = uColCountA; n < uNewColCount; ++n)\r
+ {\r
+ for (unsigned uSeqIndexA = 0; uSeqIndexA < uSeqCountA; ++uSeqIndexA)\r
+ msaCombined.SetChar(uSeqIndexA, uColIndexCombined + n, '.');\r
+ }\r
+\r
+ for (unsigned n = 0; n < uColCountB; ++n)\r
+ {\r
+ for (unsigned uSeqIndexB = 0; uSeqIndexB < uSeqCountB; ++uSeqIndexB)\r
+ {\r
+ char c = msaB.GetChar(uSeqIndexB, uColIndexB + n);\r
+ c = UnalignChar(c);\r
+ msaCombined.SetChar(uSeqCountA + uSeqIndexB, uColIndexCombined + n, c);\r
+ }\r
+ }\r
+ for (unsigned n = uColCountB; n < uNewColCount; ++n)\r
+ {\r
+ for (unsigned uSeqIndexB = 0; uSeqIndexB < uSeqCountB; ++uSeqIndexB)\r
+ msaCombined.SetChar(uSeqCountA + uSeqIndexB, uColIndexCombined + n, '.');\r
+ }\r
+\r
+ uColIndexCombined += uNewColCount;\r
+ uColIndexA += uColCountA;\r
+ uColIndexB += uColCountB;\r
+ }\r
+\r
+static void AppendMatch(const MSA &msaA, unsigned &uColIndexA, const MSA &msaB,\r
+ unsigned &uColIndexB, unsigned uSeqCountA, unsigned uSeqCountB,\r
+ MSA &msaCombined, unsigned &uColIndexCombined)\r
+ {\r
+#if TRACE\r
+ Log("AppendMatch ColIxA=%u ColIxB=%u ColIxCmb=%u\n",\r
+ uColIndexA, uColIndexB, uColIndexCombined);\r
+#endif\r
+\r
+ for (unsigned uSeqIndexA = 0; uSeqIndexA < uSeqCountA; ++uSeqIndexA)\r
+ {\r
+ char c = msaA.GetChar(uSeqIndexA, uColIndexA);\r
+ msaCombined.SetChar(uSeqIndexA, uColIndexCombined, c);\r
+ }\r
+\r
+ for (unsigned uSeqIndexB = 0; uSeqIndexB < uSeqCountB; ++uSeqIndexB)\r
+ {\r
+ char c = msaB.GetChar(uSeqIndexB, uColIndexB);\r
+ msaCombined.SetChar(uSeqCountA + uSeqIndexB, uColIndexCombined, c);\r
+ }\r
+\r
+ ++uColIndexA;\r
+ ++uColIndexB;\r
+ ++uColIndexCombined;\r
+ }\r
+\r
+void AlignTwoMSAsGivenPath(const PWPath &Path, const MSA &msaA, const MSA &msaB,\r
+ MSA &msaCombined)\r
+ {\r
+ msaCombined.Clear();\r
+\r
+#if TRACE\r
+ Log("FastAlignProfiles\n");\r
+ Log("Template A:\n");\r
+ msaA.LogMe();\r
+ Log("Template B:\n");\r
+ msaB.LogMe();\r
+#endif\r
+\r
+ const unsigned uColCountA = msaA.GetColCount();\r
+ const unsigned uColCountB = msaB.GetColCount();\r
+\r
+ const unsigned uSeqCountA = msaA.GetSeqCount();\r
+ const unsigned uSeqCountB = msaB.GetSeqCount();\r
+\r
+ msaCombined.SetSeqCount(uSeqCountA + uSeqCountB);\r
+\r
+// Copy sequence names into combined MSA\r
+ for (unsigned uSeqIndexA = 0; uSeqIndexA < uSeqCountA; ++uSeqIndexA)\r
+ {\r
+ msaCombined.SetSeqName(uSeqIndexA, msaA.GetSeqName(uSeqIndexA));\r
+ msaCombined.SetSeqId(uSeqIndexA, msaA.GetSeqId(uSeqIndexA));\r
+ }\r
+\r
+ for (unsigned uSeqIndexB = 0; uSeqIndexB < uSeqCountB; ++uSeqIndexB)\r
+ {\r
+ msaCombined.SetSeqName(uSeqCountA + uSeqIndexB, msaB.GetSeqName(uSeqIndexB));\r
+ msaCombined.SetSeqId(uSeqCountA + uSeqIndexB, msaB.GetSeqId(uSeqIndexB));\r
+ }\r
+\r
+ unsigned uColIndexA = 0;\r
+ unsigned uColIndexB = 0;\r
+ unsigned uColIndexCombined = 0;\r
+ const unsigned uEdgeCount = Path.GetEdgeCount();\r
+ for (unsigned uEdgeIndex = 0; uEdgeIndex < uEdgeCount; ++uEdgeIndex)\r
+ {\r
+ const PWEdge &Edge = Path.GetEdge(uEdgeIndex);\r
+#if TRACE\r
+ Log("\nEdge %u %c%u.%u\n",\r
+ uEdgeIndex,\r
+ Edge.cType,\r
+ Edge.uPrefixLengthA,\r
+ Edge.uPrefixLengthB);\r
+#endif\r
+ const char cType = Edge.cType;\r
+ const unsigned uPrefixLengthA = Edge.uPrefixLengthA;\r
+ unsigned uColCountA = 0;\r
+ if (uPrefixLengthA > 0)\r
+ {\r
+ const unsigned uNodeIndexA = uPrefixLengthA - 1;\r
+ const unsigned uTplColIndexA = uNodeIndexA;\r
+ if (uTplColIndexA > uColIndexA)\r
+ uColCountA = uTplColIndexA - uColIndexA;\r
+ }\r
+\r
+ const unsigned uPrefixLengthB = Edge.uPrefixLengthB;\r
+ unsigned uColCountB = 0;\r
+ if (uPrefixLengthB > 0)\r
+ {\r
+ const unsigned uNodeIndexB = uPrefixLengthB - 1;\r
+ const unsigned uTplColIndexB = uNodeIndexB;\r
+ if (uTplColIndexB > uColIndexB)\r
+ uColCountB = uTplColIndexB - uColIndexB;\r
+ }\r
+\r
+// TODO: This code looks like a hangover from HMM estimation -- can we delete it?\r
+ assert(uColCountA == 0);\r
+ assert(uColCountB == 0);\r
+ AppendTplInserts(msaA, uColIndexA, uColCountA, msaB, uColIndexB, uColCountB,\r
+ uSeqCountA, uSeqCountB, msaCombined, uColIndexCombined);\r
+\r
+ switch (cType)\r
+ {\r
+ case 'M':\r
+ {\r
+ assert(uPrefixLengthA > 0);\r
+ assert(uPrefixLengthB > 0);\r
+ const unsigned uColA = uPrefixLengthA - 1;\r
+ const unsigned uColB = uPrefixLengthB - 1;\r
+ assert(uColIndexA == uColA);\r
+ assert(uColIndexB == uColB);\r
+ AppendMatch(msaA, uColIndexA, msaB, uColIndexB, uSeqCountA, uSeqCountB,\r
+ msaCombined, uColIndexCombined);\r
+ break;\r
+ }\r
+ case 'D':\r
+ {\r
+ assert(uPrefixLengthA > 0);\r
+ const unsigned uColA = uPrefixLengthA - 1;\r
+ assert(uColIndexA == uColA);\r
+ AppendDelete(msaA, uColIndexA, uSeqCountA, uSeqCountB, msaCombined, uColIndexCombined);\r
+ break;\r
+ }\r
+ case 'I':\r
+ {\r
+ assert(uPrefixLengthB > 0);\r
+ const unsigned uColB = uPrefixLengthB - 1;\r
+ assert(uColIndexB == uColB);\r
+ AppendInsert(msaB, uColIndexB, uSeqCountA, uSeqCountB, msaCombined, uColIndexCombined);\r
+ break;\r
+ }\r
+ default:\r
+ assert(false);\r
+ }\r
+ }\r
+ unsigned uInsertColCountA = uColCountA - uColIndexA;\r
+ unsigned uInsertColCountB = uColCountB - uColIndexB;\r
+\r
+// TODO: This code looks like a hangover from HMM estimation -- can we delete it?\r
+ assert(uInsertColCountA == 0);\r
+ assert(uInsertColCountB == 0);\r
+ AppendTplInserts(msaA, uColIndexA, uInsertColCountA, msaB, uColIndexB,\r
+ uInsertColCountB, uSeqCountA, uSeqCountB, msaCombined, uColIndexCombined);\r
+\r
+ assert(msaCombined.GetColCount() == uEdgeCount);\r
+ }\r
+\r
+static const ProfPos PPStart =\r
+ {\r
+ false, //m_bAllGaps;\r
+ { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // m_uSortOrder[21];\r
+ { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // m_fcCounts[20];\r
+ 1.0, // m_LL;\r
+ 0.0, // m_LG;\r
+ 0.0, // m_GL;\r
+ 0.0, // m_GG;\r
+ { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // m_ALScores\r
+ 0, // m_uResidueGroup;\r
+ 1.0, // m_fOcc;\r
+ 0.0, // m_fcStartOcc;\r
+ 0.0, // m_fcEndOcc;\r
+ 0.0, // m_scoreGapOpen;\r
+ 0.0, // m_scoreGapClose;\r
+ };\r
+\r
+// MM\r
+// Ai\961 Ai Out\r
+// X X LL LL\r
+// X - LG LG\r
+// - X GL GL\r
+// - - GG GG\r
+// \r
+// Bj\961 Bj\r
+// X X LL LL\r
+// X - LG LG\r
+// - X GL GL\r
+// - - GG GG\r
+static void SetGapsMM(\r
+ const ProfPos *PA, unsigned uPrefixLengthA, WEIGHT wA,\r
+ const ProfPos *PB, unsigned uPrefixLengthB, WEIGHT wB,\r
+ ProfPos *POut, unsigned uColIndexOut)\r
+ {\r
+ const ProfPos &PPA = uPrefixLengthA > 0 ? PA[uPrefixLengthA-1] : PPStart;\r
+ const ProfPos &PPB = uPrefixLengthB > 0 ? PB[uPrefixLengthB-1] : PPStart;\r
+ ProfPos &PPO = POut[uColIndexOut];\r
+\r
+ PPO.m_LL = wA*PPA.m_LL + wB*PPB.m_LL;\r
+ PPO.m_LG = wA*PPA.m_LG + wB*PPB.m_LG;\r
+ PPO.m_GL = wA*PPA.m_GL + wB*PPB.m_GL;\r
+ PPO.m_GG = wA*PPA.m_GG + wB*PPB.m_GG;\r
+ }\r
+\r
+// MD\r
+// Ai\961 Ai Out\r
+// X X LL LL\r
+// X - LG LG\r
+// - X GL GL\r
+// - - GG GG\r
+// \r
+// Bj (-)\r
+// X - ?L LG\r
+// - - ?G GG\r
+static void SetGapsMD(\r
+ const ProfPos *PA, unsigned uPrefixLengthA, WEIGHT wA,\r
+ const ProfPos *PB, unsigned uPrefixLengthB, WEIGHT wB,\r
+ ProfPos *POut, unsigned uColIndexOut)\r
+ {\r
+ const ProfPos &PPA = uPrefixLengthA > 0 ? PA[uPrefixLengthA-1] : PPStart;\r
+ const ProfPos &PPB = uPrefixLengthB > 0 ? PB[uPrefixLengthB-1] : PPStart;\r
+ ProfPos &PPO = POut[uColIndexOut];\r
+\r
+ PPO.m_LL = wA*PPA.m_LL;\r
+ PPO.m_LG = wA*PPA.m_LG + wB*(PPB.m_LL + PPB.m_GL);\r
+ PPO.m_GL = wA*PPA.m_GL;\r
+ PPO.m_GG = wA*PPA.m_GG + wB*(PPB.m_LG + PPB.m_GG);\r
+ }\r
+\r
+// DD\r
+// Ai\961 Ai Out\r
+// X X LL LL\r
+// X - LG LG\r
+// - X GL GL\r
+// - - GG GG\r
+// \r
+// (-) (-)\r
+// - - ?? GG\r
+static void SetGapsDD(\r
+ const ProfPos *PA, unsigned uPrefixLengthA, WEIGHT wA,\r
+ const ProfPos *PB, unsigned uPrefixLengthB, WEIGHT wB,\r
+ ProfPos *POut, unsigned uColIndexOut)\r
+ {\r
+ const ProfPos &PPA = uPrefixLengthA > 0 ? PA[uPrefixLengthA-1] : PPStart;\r
+ ProfPos &PPO = POut[uColIndexOut];\r
+\r
+ PPO.m_LL = wA*PPA.m_LL;\r
+ PPO.m_LG = wA*PPA.m_LG;\r
+ PPO.m_GL = wA*PPA.m_GL;\r
+ PPO.m_GG = wA*PPA.m_GG + wB;\r
+ }\r
+\r
+// MI\r
+// Ai (-) Out\r
+// X - ?L LG\r
+// - - ?G GG\r
+\r
+// Bj\961 Bj\r
+// X X LL LL\r
+// X - LG LG\r
+// - X GL GL\r
+// - - GG GG\r
+static void SetGapsMI(\r
+ const ProfPos *PA, unsigned uPrefixLengthA, WEIGHT wA,\r
+ const ProfPos *PB, unsigned uPrefixLengthB, WEIGHT wB,\r
+ ProfPos *POut, unsigned uColIndexOut)\r
+ {\r
+ const ProfPos &PPA = uPrefixLengthA > 0 ? PA[uPrefixLengthA-1] : PPStart;\r
+ const ProfPos &PPB = uPrefixLengthB > 0 ? PB[uPrefixLengthB-1] : PPStart;\r
+ ProfPos &PPO = POut[uColIndexOut];\r
+\r
+ PPO.m_LL = wB*PPB.m_LL;\r
+ PPO.m_LG = wB*PPB.m_LG + wA*(PPA.m_LL + PPA.m_GL);\r
+ PPO.m_GL = wB*PPB.m_GL;\r
+ PPO.m_GG = wB*PPB.m_GG + wA*(PPA.m_LG + PPA.m_GG);\r
+ }\r
+\r
+// DM\r
+// Ai\961 Ai Out\r
+// X X LL LL\r
+// X - LG LG\r
+// - X GL GL\r
+// - - GG GG\r
+// \r
+// (-) Bj \r
+// - X ?L GL\r
+// - - ?G GG\r
+static void SetGapsDM(\r
+ const ProfPos *PA, unsigned uPrefixLengthA, WEIGHT wA,\r
+ const ProfPos *PB, unsigned uPrefixLengthB, WEIGHT wB,\r
+ ProfPos *POut, unsigned uColIndexOut)\r
+ {\r
+ const ProfPos &PPA = uPrefixLengthA > 0 ? PA[uPrefixLengthA-1] : PPStart;\r
+ const ProfPos &PPB = uPrefixLengthB > 0 ? PB[uPrefixLengthB-1] : PPStart;\r
+ ProfPos &PPO = POut[uColIndexOut];\r
+\r
+ PPO.m_LL = wA*PPA.m_LL;\r
+ PPO.m_LG = wA*PPA.m_LG;\r
+ PPO.m_GL = wA*PPA.m_GL + wB*(PPB.m_LL + PPB.m_GL);\r
+ PPO.m_GG = wA*PPA.m_GG + wB*(PPB.m_LG + PPB.m_GG);\r
+ }\r
+\r
+// IM\r
+// (-) Ai Out \r
+// - X ?L GL\r
+// - - ?G GG\r
+\r
+// Bj\961 Bj\r
+// X X LL LL\r
+// X - LG LG\r
+// - X GL GL\r
+// - - GG GG\r
+static void SetGapsIM(\r
+ const ProfPos *PA, unsigned uPrefixLengthA, WEIGHT wA,\r
+ const ProfPos *PB, unsigned uPrefixLengthB, WEIGHT wB,\r
+ ProfPos *POut, unsigned uColIndexOut)\r
+ {\r
+ const ProfPos &PPA = uPrefixLengthA > 0 ? PA[uPrefixLengthA-1] : PPStart;\r
+ const ProfPos &PPB = uPrefixLengthB > 0 ? PB[uPrefixLengthB-1] : PPStart;\r
+ ProfPos &PPO = POut[uColIndexOut];\r
+\r
+ PPO.m_LL = wB*PPB.m_LL;\r
+ PPO.m_LG = wB*PPB.m_LG;\r
+ PPO.m_GL = wB*PPB.m_GL + wA*(PPA.m_LL + PPA.m_GL);\r
+ PPO.m_GG = wB*PPB.m_GG + wA*(PPA.m_LG + PPA.m_GG);\r
+ }\r
+\r
+// ID\r
+// (-) Ai Out\r
+// - X ?L GL\r
+// - - ?G GG\r
+\r
+// Bj (-)\r
+// X - ?L LG\r
+// - - ?G GG\r
+static void SetGapsID(\r
+ const ProfPos *PA, unsigned uPrefixLengthA, WEIGHT wA,\r
+ const ProfPos *PB, unsigned uPrefixLengthB, WEIGHT wB,\r
+ ProfPos *POut, unsigned uColIndexOut)\r
+ {\r
+ const ProfPos &PPA = uPrefixLengthA > 0 ? PA[uPrefixLengthA-1] : PPStart;\r
+ const ProfPos &PPB = uPrefixLengthB > 0 ? PB[uPrefixLengthB-1] : PPStart;\r
+ ProfPos &PPO = POut[uColIndexOut];\r
+\r
+ PPO.m_LL = 0;\r
+ PPO.m_LG = wB*PPB.m_GL + wB*PPB.m_LL;\r
+ PPO.m_GL = wA*PPA.m_GL + wA*PPA.m_LL;\r
+ PPO.m_GG = wA*(PPA.m_LG + PPA.m_GG) + wB*(PPB.m_LG + PPB.m_GG);\r
+ }\r
+\r
+// DI\r
+// Ai (-) Out\r
+// X - ?L LG\r
+// - - ?G GG\r
+\r
+// (-) Bj\r
+// - X ?L GL\r
+// - - ?G GG\r
+static void SetGapsDI(\r
+ const ProfPos *PA, unsigned uPrefixLengthA, WEIGHT wA,\r
+ const ProfPos *PB, unsigned uPrefixLengthB, WEIGHT wB,\r
+ ProfPos *POut, unsigned uColIndexOut)\r
+ {\r
+ const ProfPos &PPA = uPrefixLengthA > 0 ? PA[uPrefixLengthA-1] : PPStart;\r
+ const ProfPos &PPB = uPrefixLengthB > 0 ? PB[uPrefixLengthB-1] : PPStart;\r
+ ProfPos &PPO = POut[uColIndexOut];\r
+\r
+ PPO.m_LL = 0;\r
+ PPO.m_LG = wA*PPA.m_GL + wA*PPA.m_LL;\r
+ PPO.m_GL = wB*PPB.m_GL + wB*PPB.m_LL;\r
+ PPO.m_GG = wA*(PPA.m_LG + PPA.m_GG) + wB*(PPB.m_LG + PPB.m_GG);\r
+ }\r
+\r
+// II\r
+// (-) (-) Out\r
+// - - ?? GG\r
+\r
+// Bj\961 Bj\r
+// X X LL LL\r
+// X - LG LG\r
+// - X GL GL\r
+// - - GG GG\r
+static void SetGapsII(\r
+ const ProfPos *PA, unsigned uPrefixLengthA, WEIGHT wA,\r
+ const ProfPos *PB, unsigned uPrefixLengthB, WEIGHT wB,\r
+ ProfPos *POut, unsigned uColIndexOut)\r
+ {\r
+ const ProfPos &PPB = uPrefixLengthB > 0 ? PB[uPrefixLengthB-1] : PPStart;\r
+ ProfPos &PPO = POut[uColIndexOut];\r
+\r
+ PPO.m_LL = wB*PPB.m_LL;\r
+ PPO.m_LG = wB*PPB.m_LG;\r
+ PPO.m_GL = wB*PPB.m_GL;\r
+ PPO.m_GG = wB*PPB.m_GG + wA;\r
+ }\r
+\r
+static void SetFreqs(\r
+ const ProfPos *PA, unsigned uPrefixLengthA, WEIGHT wA,\r
+ const ProfPos *PB, unsigned uPrefixLengthB, WEIGHT wB,\r
+ ProfPos *POut, unsigned uColIndexOut)\r
+ {\r
+ const ProfPos &PPA = uPrefixLengthA > 0 ? PA[uPrefixLengthA-1] : PPStart;\r
+ const ProfPos &PPB = uPrefixLengthB > 0 ? PB[uPrefixLengthB-1] : PPStart;\r
+ ProfPos &PPO = POut[uColIndexOut];\r
+\r
+ if (g_bNormalizeCounts)\r
+ {\r
+ const FCOUNT fA = PPA.m_fOcc*wA/(wA + wB);\r
+ const FCOUNT fB = PPB.m_fOcc*wB/(wA + wB);\r
+ FCOUNT fTotal = 0;\r
+ for (unsigned i = 0; i < 20; ++i)\r
+ {\r
+ const FCOUNT f = fA*PPA.m_fcCounts[i] + fB*PPB.m_fcCounts[i];\r
+ PPO.m_fcCounts[i] = f;\r
+ fTotal += f;\r
+ }\r
+ if (fTotal > 0)\r
+ for (unsigned i = 0; i < 20; ++i)\r
+ PPO.m_fcCounts[i] /= fTotal;\r
+ }\r
+ else\r
+ {\r
+ for (unsigned i = 0; i < 20; ++i)\r
+ PPO.m_fcCounts[i] = wA*PPA.m_fcCounts[i] + wB*PPB.m_fcCounts[i];\r
+ }\r
+ }\r
+\r
+void AlignTwoProfsGivenPath(const PWPath &Path,\r
+ const ProfPos *PA, unsigned uPrefixLengthA, WEIGHT wA,\r
+ const ProfPos *PB, unsigned uPrefixLengthB, WEIGHT wB,\r
+ ProfPos **ptrPOut, unsigned *ptruLengthOut)\r
+ {\r
+#if TRACE\r
+ Log("AlignTwoProfsGivenPath wA=%.3g wB=%.3g Path=\n", wA, wB);\r
+ Path.LogMe();\r
+#endif\r
+ assert(BTEq(wA + wB, 1.0));\r
+\r
+ unsigned uColIndexA = 0;\r
+ unsigned uColIndexB = 0;\r
+ unsigned uColIndexOut = 0;\r
+ const unsigned uEdgeCount = Path.GetEdgeCount();\r
+ ProfPos *POut = new ProfPos[uEdgeCount];\r
+ char cPrevType = 'M';\r
+ for (unsigned uEdgeIndex = 0; uEdgeIndex < uEdgeCount; ++uEdgeIndex)\r
+ {\r
+ const PWEdge &Edge = Path.GetEdge(uEdgeIndex);\r
+ const char cType = Edge.cType;\r
+\r
+ const unsigned uPrefixLengthA = Edge.uPrefixLengthA;\r
+ const unsigned uPrefixLengthB = Edge.uPrefixLengthB;\r
+\r
+#if TRACE\r
+ Log("\nEdge %u %c%u.%u ColA=%u ColB=%u\n",\r
+ uEdgeIndex,\r
+ Edge.cType,\r
+ Edge.uPrefixLengthA,\r
+ Edge.uPrefixLengthB,\r
+ uColIndexA,\r
+ uColIndexB);\r
+#endif\r
+\r
+ POut[uColIndexOut].m_bAllGaps = false;\r
+ switch (cType)\r
+ {\r
+ case 'M':\r
+ {\r
+ assert(uPrefixLengthA > 0);\r
+ assert(uPrefixLengthB > 0);\r
+ SetFreqs(\r
+ PA, uPrefixLengthA, wA,\r
+ PB, uPrefixLengthB, wB,\r
+ POut, uColIndexOut);\r
+ switch (cPrevType)\r
+ {\r
+ case 'M':\r
+ SetGapsMM(\r
+ PA, uPrefixLengthA, wA,\r
+ PB, uPrefixLengthB, wB,\r
+ POut, uColIndexOut);\r
+ break;\r
+ case 'D':\r
+ SetGapsDM(\r
+ PA, uPrefixLengthA, wA,\r
+ PB, uPrefixLengthB, wB,\r
+ POut, uColIndexOut);\r
+ break;\r
+ case 'I':\r
+ SetGapsIM(\r
+ PA, uPrefixLengthA, wA,\r
+ PB, uPrefixLengthB, wB,\r
+ POut, uColIndexOut);\r
+ break;\r
+ default:\r
+ Quit("Bad cPrevType");\r
+ }\r
+ ++uColIndexA;\r
+ ++uColIndexB;\r
+ ++uColIndexOut;\r
+ break;\r
+ }\r
+ case 'D':\r
+ {\r
+ assert(uPrefixLengthA > 0);\r
+ SetFreqs(\r
+ PA, uPrefixLengthA, wA,\r
+ PB, uPrefixLengthB, 0,\r
+ POut, uColIndexOut);\r
+ switch (cPrevType)\r
+ {\r
+ case 'M':\r
+ SetGapsMD(\r
+ PA, uPrefixLengthA, wA,\r
+ PB, uPrefixLengthB, wB,\r
+ POut, uColIndexOut);\r
+ break;\r
+ case 'D':\r
+ SetGapsDD(\r
+ PA, uPrefixLengthA, wA,\r
+ PB, uPrefixLengthB, wB,\r
+ POut, uColIndexOut);\r
+ break;\r
+ case 'I':\r
+ SetGapsID(\r
+ PA, uPrefixLengthA, wA,\r
+ PB, uPrefixLengthB, wB,\r
+ POut, uColIndexOut);\r
+ break;\r
+ default:\r
+ Quit("Bad cPrevType");\r
+ }\r
+ ++uColIndexA;\r
+ ++uColIndexOut;\r
+ break;\r
+ }\r
+ case 'I':\r
+ {\r
+ assert(uPrefixLengthB > 0);\r
+ SetFreqs(\r
+ PA, uPrefixLengthA, 0,\r
+ PB, uPrefixLengthB, wB,\r
+ POut, uColIndexOut);\r
+ switch (cPrevType)\r
+ {\r
+ case 'M':\r
+ SetGapsMI(\r
+ PA, uPrefixLengthA, wA,\r
+ PB, uPrefixLengthB, wB,\r
+ POut, uColIndexOut);\r
+ break;\r
+ case 'D':\r
+ SetGapsDI(\r
+ PA, uPrefixLengthA, wA,\r
+ PB, uPrefixLengthB, wB,\r
+ POut, uColIndexOut);\r
+ break;\r
+ case 'I':\r
+ SetGapsII(\r
+ PA, uPrefixLengthA, wA,\r
+ PB, uPrefixLengthB, wB,\r
+ POut, uColIndexOut);\r
+ break;\r
+ default:\r
+ Quit("Bad cPrevType");\r
+ }\r
+ ++uColIndexB;\r
+ ++uColIndexOut;\r
+ break;\r
+ }\r
+ default:\r
+ assert(false);\r
+ }\r
+ cPrevType = cType;\r
+ }\r
+ assert(uColIndexOut == uEdgeCount);\r
+\r
+ ProfScoresFromFreqs(POut, uEdgeCount);\r
+ ValidateProf(POut, uEdgeCount);\r
+\r
+ *ptrPOut = POut;\r
+ *ptruLengthOut = uEdgeCount;\r
+\r
+#if TRACE\r
+ Log("AlignTwoProfsGivenPath:\n");\r
+ ListProfile(POut, uEdgeCount, 0);\r
+#endif\r
+ }\r
git://source.jalview.org / jabaws.git / blobdiff