--- /dev/null
+#include "muscle.h"\r
+#include "msa.h"\r
+#include "profile.h"\r
+#include "objscore.h"\r
+\r
+#define TRACE 0\r
+#define TRACE_SEQPAIR 0\r
+#define TEST_SPFAST 0\r
+\r
+extern SCOREMATRIX VTML_LA;\r
+extern SCOREMATRIX PAM200;\r
+extern SCOREMATRIX PAM200NoCenter;\r
+extern SCOREMATRIX VTML_SP;\r
+extern SCOREMATRIX VTML_SPNoCenter;\r
+extern SCOREMATRIX NUC_SP;\r
+\r
+SCORE g_SPScoreLetters;\r
+SCORE g_SPScoreGaps;\r
+\r
+static SCORE TermGapScore(bool Gap)\r
+ {\r
+ switch (g_TermGaps)\r
+ {\r
+ case TERMGAPS_Full:\r
+ return 0;\r
+\r
+ case TERMGAPS_Half:\r
+ if (Gap)\r
+ return g_scoreGapOpen/2;\r
+ return 0;\r
+\r
+ case TERMGAPS_Ext:\r
+ if (Gap)\r
+ return g_scoreGapExtend;\r
+ return 0;\r
+ }\r
+ Quit("TermGapScore?!");\r
+ return 0;\r
+ }\r
+\r
+SCORE ScoreSeqPairLetters(const MSA &msa1, unsigned uSeqIndex1,\r
+ const MSA &msa2, unsigned uSeqIndex2)\r
+ {\r
+ const unsigned uColCount = msa1.GetColCount();\r
+ const unsigned uColCount2 = msa2.GetColCount();\r
+ if (uColCount != uColCount2)\r
+ Quit("ScoreSeqPairLetters, different lengths");\r
+\r
+#if TRACE_SEQPAIR\r
+ {\r
+ Log("\n");\r
+ Log("ScoreSeqPairLetters\n");\r
+ MSA msaTmp;\r
+ msaTmp.SetSize(2, uColCount);\r
+ msaTmp.CopySeq(0, msa1, uSeqIndex1);\r
+ msaTmp.CopySeq(1, msa2, uSeqIndex2);\r
+ msaTmp.LogMe();\r
+ }\r
+#endif\r
+\r
+ SCORE scoreLetters = 0;\r
+ SCORE scoreGaps = 0;\r
+ bool bGapping1 = false;\r
+ bool bGapping2 = false;\r
+\r
+ unsigned uColStart = 0;\r
+ bool bLeftTermGap = false;\r
+ for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)\r
+ {\r
+ bool bGap1 = msa1.IsGap(uSeqIndex1, uColIndex);\r
+ bool bGap2 = msa2.IsGap(uSeqIndex2, uColIndex);\r
+ if (!bGap1 || !bGap2)\r
+ {\r
+ if (bGap1 || bGap2)\r
+ bLeftTermGap = true;\r
+ uColStart = uColIndex;\r
+ break;\r
+ }\r
+ }\r
+\r
+ unsigned uColEnd = uColCount - 1;\r
+ bool bRightTermGap = false;\r
+ for (int iColIndex = (int) uColCount - 1; iColIndex >= 0; --iColIndex)\r
+ {\r
+ bool bGap1 = msa1.IsGap(uSeqIndex1, iColIndex);\r
+ bool bGap2 = msa2.IsGap(uSeqIndex2, iColIndex);\r
+ if (!bGap1 || !bGap2)\r
+ {\r
+ if (bGap1 || bGap2)\r
+ bRightTermGap = true;\r
+ uColEnd = (unsigned) iColIndex;\r
+ break;\r
+ }\r
+ }\r
+\r
+#if TRACE_SEQPAIR\r
+ Log("LeftTermGap=%d RightTermGap=%d\n", bLeftTermGap, bRightTermGap);\r
+#endif\r
+\r
+ for (unsigned uColIndex = uColStart; uColIndex <= uColEnd; ++uColIndex)\r
+ {\r
+ unsigned uLetter1 = msa1.GetLetterEx(uSeqIndex1, uColIndex);\r
+ if (uLetter1 >= g_AlphaSize)\r
+ continue;\r
+ unsigned uLetter2 = msa2.GetLetterEx(uSeqIndex2, uColIndex);\r
+ if (uLetter2 >= g_AlphaSize)\r
+ continue;\r
+\r
+ SCORE scoreMatch = (*g_ptrScoreMatrix)[uLetter1][uLetter2];\r
+ scoreLetters += scoreMatch;\r
+ }\r
+ return scoreLetters;\r
+ }\r
+\r
+SCORE ScoreSeqPairGaps(const MSA &msa1, unsigned uSeqIndex1,\r
+ const MSA &msa2, unsigned uSeqIndex2)\r
+ {\r
+ const unsigned uColCount = msa1.GetColCount();\r
+ const unsigned uColCount2 = msa2.GetColCount();\r
+ if (uColCount != uColCount2)\r
+ Quit("ScoreSeqPairGaps, different lengths");\r
+\r
+#if TRACE_SEQPAIR\r
+ {\r
+ Log("\n");\r
+ Log("ScoreSeqPairGaps\n");\r
+ MSA msaTmp;\r
+ msaTmp.SetSize(2, uColCount);\r
+ msaTmp.CopySeq(0, msa1, uSeqIndex1);\r
+ msaTmp.CopySeq(1, msa2, uSeqIndex2);\r
+ msaTmp.LogMe();\r
+ }\r
+#endif\r
+\r
+ SCORE scoreGaps = 0;\r
+ bool bGapping1 = false;\r
+ bool bGapping2 = false;\r
+\r
+ unsigned uColStart = 0;\r
+ bool bLeftTermGap = false;\r
+ for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)\r
+ {\r
+ bool bGap1 = msa1.IsGap(uSeqIndex1, uColIndex);\r
+ bool bGap2 = msa2.IsGap(uSeqIndex2, uColIndex);\r
+ if (!bGap1 || !bGap2)\r
+ {\r
+ if (bGap1 || bGap2)\r
+ bLeftTermGap = true;\r
+ uColStart = uColIndex;\r
+ break;\r
+ }\r
+ }\r
+\r
+ unsigned uColEnd = uColCount - 1;\r
+ bool bRightTermGap = false;\r
+ for (int iColIndex = (int) uColCount - 1; iColIndex >= 0; --iColIndex)\r
+ {\r
+ bool bGap1 = msa1.IsGap(uSeqIndex1, iColIndex);\r
+ bool bGap2 = msa2.IsGap(uSeqIndex2, iColIndex);\r
+ if (!bGap1 || !bGap2)\r
+ {\r
+ if (bGap1 || bGap2)\r
+ bRightTermGap = true;\r
+ uColEnd = (unsigned) iColIndex;\r
+ break;\r
+ }\r
+ }\r
+\r
+#if TRACE_SEQPAIR\r
+ Log("LeftTermGap=%d RightTermGap=%d\n", bLeftTermGap, bRightTermGap);\r
+#endif\r
+\r
+ for (unsigned uColIndex = uColStart; uColIndex <= uColEnd; ++uColIndex)\r
+ {\r
+ bool bGap1 = msa1.IsGap(uSeqIndex1, uColIndex);\r
+ bool bGap2 = msa2.IsGap(uSeqIndex2, uColIndex);\r
+\r
+ if (bGap1 && bGap2)\r
+ continue;\r
+\r
+ if (bGap1)\r
+ {\r
+ if (!bGapping1)\r
+ {\r
+#if TRACE_SEQPAIR\r
+ Log("Gap open seq 1 col %d\n", uColIndex);\r
+#endif\r
+ if (uColIndex == uColStart)\r
+ scoreGaps += TermGapScore(true);\r
+ else\r
+ scoreGaps += g_scoreGapOpen;\r
+ bGapping1 = true;\r
+ }\r
+ else\r
+ scoreGaps += g_scoreGapExtend;\r
+ continue;\r
+ }\r
+\r
+ else if (bGap2)\r
+ {\r
+ if (!bGapping2)\r
+ {\r
+#if TRACE_SEQPAIR\r
+ Log("Gap open seq 2 col %d\n", uColIndex);\r
+#endif\r
+ if (uColIndex == uColStart)\r
+ scoreGaps += TermGapScore(true);\r
+ else\r
+ scoreGaps += g_scoreGapOpen;\r
+ bGapping2 = true;\r
+ }\r
+ else\r
+ scoreGaps += g_scoreGapExtend;\r
+ continue;\r
+ }\r
+\r
+ bGapping1 = false;\r
+ bGapping2 = false;\r
+ }\r
+\r
+ if (bGapping1 || bGapping2)\r
+ {\r
+ scoreGaps -= g_scoreGapOpen;\r
+ scoreGaps += TermGapScore(true);\r
+ }\r
+ return scoreGaps;\r
+ }\r
+\r
+// The usual sum-of-pairs objective score: sum the score\r
+// of the alignment of each pair of sequences.\r
+SCORE ObjScoreSP(const MSA &msa, SCORE MatchScore[])\r
+ {\r
+#if TRACE\r
+ Log("==================ObjScoreSP==============\n");\r
+ Log("msa=\n");\r
+ msa.LogMe();\r
+#endif\r
+ g_SPScoreLetters = 0;\r
+ g_SPScoreGaps = 0;\r
+\r
+ if (0 != MatchScore)\r
+ {\r
+ const unsigned uColCount = msa.GetColCount();\r
+ for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)\r
+ MatchScore[uColIndex] = 0;\r
+ }\r
+\r
+ const unsigned uSeqCount = msa.GetSeqCount();\r
+ SCORE scoreTotal = 0;\r
+ unsigned uPairCount = 0;\r
+#if TRACE\r
+ Log("Seq1 Seq2 wt1 wt2 Letters Gaps Unwt.Score Wt.Score Total\n");\r
+ Log("---- ---- ------ ------ ---------- ---------- ---------- ---------- ----------\n");\r
+#endif\r
+ for (unsigned uSeqIndex1 = 0; uSeqIndex1 < uSeqCount; ++uSeqIndex1)\r
+ {\r
+ const WEIGHT w1 = msa.GetSeqWeight(uSeqIndex1);\r
+ for (unsigned uSeqIndex2 = uSeqIndex1 + 1; uSeqIndex2 < uSeqCount; ++uSeqIndex2)\r
+ {\r
+ const WEIGHT w2 = msa.GetSeqWeight(uSeqIndex2);\r
+ const WEIGHT w = w1*w2;\r
+\r
+ SCORE scoreLetters = ScoreSeqPairLetters(msa, uSeqIndex1, msa, uSeqIndex2);\r
+ SCORE scoreGaps = ScoreSeqPairGaps(msa, uSeqIndex1, msa, uSeqIndex2);\r
+ SCORE scorePair = scoreLetters + scoreGaps;\r
+ ++uPairCount;\r
+\r
+ scoreTotal += w*scorePair;\r
+\r
+ g_SPScoreLetters += w*scoreLetters;\r
+ g_SPScoreGaps += w*scoreGaps;\r
+#if TRACE\r
+ Log("%4d %4d %6.3f %6.3f %10.2f %10.2f %10.2f %10.2f %10.2f >%s >%s\n",\r
+ uSeqIndex1,\r
+ uSeqIndex2,\r
+ w1,\r
+ w2,\r
+ scoreLetters,\r
+ scoreGaps,\r
+ scorePair,\r
+ scorePair*w1*w2,\r
+ scoreTotal,\r
+ msa.GetSeqName(uSeqIndex1),\r
+ msa.GetSeqName(uSeqIndex2));\r
+#endif\r
+ }\r
+ }\r
+#if TEST_SPFAST\r
+ {\r
+ SCORE f = ObjScoreSPFast(msa);\r
+ Log("Fast = %.6g\n", f);\r
+ Log("Brute = %.6g\n", scoreTotal);\r
+ if (BTEq(f, scoreTotal))\r
+ Log("Agree\n");\r
+ else\r
+ Log("** DISAGREE **\n");\r
+ }\r
+#endif\r
+// return scoreTotal / uPairCount;\r
+ return scoreTotal;\r
+ }\r
+\r
+// Objective score defined as the dynamic programming score.\r
+// Input is two alignments, which must be of the same length.\r
+// Result is the same profile-profile score that is optimized\r
+// by dynamic programming.\r
+SCORE ObjScoreDP(const MSA &msa1, const MSA &msa2, SCORE MatchScore[])\r
+ {\r
+ const unsigned uColCount = msa1.GetColCount();\r
+ if (msa2.GetColCount() != uColCount)\r
+ Quit("ObjScoreDP, must be same length");\r
+\r
+ const unsigned uColCount1 = msa1.GetColCount();\r
+ const unsigned uColCount2 = msa2.GetColCount();\r
+\r
+ const ProfPos *PA = ProfileFromMSA(msa1);\r
+ const ProfPos *PB = ProfileFromMSA(msa2);\r
+\r
+ return ObjScoreDP_Profs(PA, PB, uColCount1, MatchScore);\r
+ }\r
+\r
+SCORE ObjScoreDP_Profs(const ProfPos *PA, const ProfPos *PB, unsigned uColCount,\r
+ SCORE MatchScore[])\r
+ {\r
+//#if TRACE\r
+// Log("Profile 1:\n");\r
+// ListProfile(PA, uColCount, &msa1);\r
+//\r
+// Log("Profile 2:\n");\r
+// ListProfile(PB, uColCount, &msa2);\r
+//#endif\r
+\r
+ SCORE scoreTotal = 0;\r
+\r
+ for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)\r
+ {\r
+ const ProfPos &PPA = PA[uColIndex];\r
+ const ProfPos &PPB = PB[uColIndex];\r
+\r
+ SCORE scoreGap = 0;\r
+ SCORE scoreMatch = 0;\r
+ // If gapped column...\r
+ if (PPA.m_bAllGaps && PPB.m_bAllGaps)\r
+ scoreGap = 0;\r
+ else if (PPA.m_bAllGaps)\r
+ {\r
+ if (uColCount - 1 == uColIndex || !PA[uColIndex+1].m_bAllGaps)\r
+ scoreGap = PPB.m_scoreGapClose;\r
+ if (0 == uColIndex || !PA[uColIndex-1].m_bAllGaps)\r
+ scoreGap += PPB.m_scoreGapOpen;\r
+ //if (0 == scoreGap)\r
+ // scoreGap = PPB.m_scoreGapExtend;\r
+ }\r
+ else if (PPB.m_bAllGaps)\r
+ {\r
+ if (uColCount - 1 == uColIndex || !PB[uColIndex+1].m_bAllGaps)\r
+ scoreGap = PPA.m_scoreGapClose;\r
+ if (0 == uColIndex || !PB[uColIndex-1].m_bAllGaps)\r
+ scoreGap += PPA.m_scoreGapOpen;\r
+ //if (0 == scoreGap)\r
+ // scoreGap = PPA.m_scoreGapExtend;\r
+ }\r
+ else\r
+ scoreMatch = ScoreProfPos2(PPA, PPB);\r
+\r
+ if (0 != MatchScore)\r
+ MatchScore[uColIndex] = scoreMatch;\r
+\r
+ scoreTotal += scoreMatch + scoreGap;\r
+\r
+ extern bool g_bTracePPScore;\r
+ extern MSA *g_ptrPPScoreMSA1;\r
+ extern MSA *g_ptrPPScoreMSA2;\r
+ if (g_bTracePPScore)\r
+ {\r
+ const MSA &msa1 = *g_ptrPPScoreMSA1;\r
+ const MSA &msa2 = *g_ptrPPScoreMSA2;\r
+ const unsigned uSeqCount1 = msa1.GetSeqCount();\r
+ const unsigned uSeqCount2 = msa2.GetSeqCount();\r
+\r
+ for (unsigned n = 0; n < uSeqCount1; ++n)\r
+ Log("%c", msa1.GetChar(n, uColIndex));\r
+ Log(" ");\r
+ for (unsigned n = 0; n < uSeqCount2; ++n)\r
+ Log("%c", msa2.GetChar(n, uColIndex));\r
+ Log(" %10.3f", scoreMatch);\r
+ if (scoreGap != 0)\r
+ Log(" %10.3f", scoreGap);\r
+ Log("\n");\r
+ }\r
+ }\r
+\r
+ delete[] PA;\r
+ delete[] PB;\r
+\r
+ return scoreTotal;\r
+ }\r
+\r
+// Objective score defined as the sum of profile-sequence\r
+// scores for each sequence in the alignment. The profile\r
+// is computed from the entire alignment, so this includes\r
+// the score of each sequence against itself. This is to\r
+// avoid recomputing the profile each time, so we reduce\r
+// complexity but introduce a questionable approximation.\r
+// The goal is to see if we can exploit the apparent\r
+// improvement in performance of log-expectation score\r
+// over the usual sum-of-pairs by optimizing this\r
+// objective score in the iterative refinement stage.\r
+SCORE ObjScorePS(const MSA &msa, SCORE MatchScore[])\r
+ {\r
+ if (g_PPScore != PPSCORE_LE)\r
+ Quit("FastScoreMSA_LASimple: LA");\r
+\r
+ const unsigned uSeqCount = msa.GetSeqCount();\r
+ const unsigned uColCount = msa.GetColCount();\r
+\r
+ const ProfPos *Prof = ProfileFromMSA(msa);\r
+\r
+ if (0 != MatchScore)\r
+ for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)\r
+ MatchScore[uColIndex] = 0;\r
+\r
+ SCORE scoreTotal = 0;\r
+ for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)\r
+ {\r
+ const WEIGHT weightSeq = msa.GetSeqWeight(uSeqIndex);\r
+ SCORE scoreSeq = 0;\r
+ for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)\r
+ {\r
+ const ProfPos &PP = Prof[uColIndex];\r
+ if (msa.IsGap(uSeqIndex, uColIndex))\r
+ {\r
+ bool bOpen = (0 == uColIndex ||\r
+ !msa.IsGap(uSeqIndex, uColIndex - 1));\r
+ bool bClose = (uColCount - 1 == uColIndex ||\r
+ !msa.IsGap(uSeqIndex, uColIndex + 1));\r
+\r
+ if (bOpen)\r
+ scoreSeq += PP.m_scoreGapOpen;\r
+ if (bClose)\r
+ scoreSeq += PP.m_scoreGapClose;\r
+ //if (!bOpen && !bClose)\r
+ // scoreSeq += PP.m_scoreGapExtend;\r
+ }\r
+ else if (msa.IsWildcard(uSeqIndex, uColIndex))\r
+ continue;\r
+ else\r
+ {\r
+ unsigned uLetter = msa.GetLetter(uSeqIndex, uColIndex);\r
+ const SCORE scoreMatch = PP.m_AAScores[uLetter];\r
+ if (0 != MatchScore)\r
+ MatchScore[uColIndex] += weightSeq*scoreMatch;\r
+ scoreSeq += scoreMatch;\r
+ }\r
+ }\r
+ scoreTotal += weightSeq*scoreSeq;\r
+ }\r
+\r
+ delete[] Prof;\r
+ return scoreTotal;\r
+ }\r
+\r
+// The XP score is the sum of the score of each pair of\r
+// sequences between two profiles which are aligned to\r
+// each other. Notice that for two given profiles aligned\r
+// in different ways, the difference in XP score must be\r
+// the same as the difference in SP score because the\r
+// score of a pair of sequences in one profile doesn't\r
+// depend on the alignment.\r
+SCORE ObjScoreXP(const MSA &msa1, const MSA &msa2)\r
+ {\r
+ const unsigned uColCount1 = msa1.GetColCount();\r
+ const unsigned uColCount2 = msa2.GetColCount();\r
+ if (uColCount1 != uColCount2)\r
+ Quit("ObjScoreXP, alignment lengths differ %u %u", uColCount1, uColCount2);\r
+\r
+ const unsigned uSeqCount1 = msa1.GetSeqCount();\r
+ const unsigned uSeqCount2 = msa2.GetSeqCount();\r
+\r
+#if TRACE\r
+ Log(" Score Weight Weight Total\n");\r
+ Log("---------- ------ ------ ----------\n");\r
+#endif\r
+\r
+ SCORE scoreTotal = 0;\r
+ unsigned uPairCount = 0;\r
+ for (unsigned uSeqIndex1 = 0; uSeqIndex1 < uSeqCount1; ++uSeqIndex1)\r
+ {\r
+ const WEIGHT w1 = msa1.GetSeqWeight(uSeqIndex1);\r
+ for (unsigned uSeqIndex2 = 0; uSeqIndex2 < uSeqCount2; ++uSeqIndex2)\r
+ {\r
+ const WEIGHT w2 = msa2.GetSeqWeight(uSeqIndex2);\r
+ const WEIGHT w = w1*w2;\r
+ SCORE scoreLetters = ScoreSeqPairLetters(msa1, uSeqIndex1, msa2, uSeqIndex2);\r
+ SCORE scoreGaps = ScoreSeqPairGaps(msa1, uSeqIndex1, msa2, uSeqIndex2);\r
+ SCORE scorePair = scoreLetters + scoreGaps;\r
+ scoreTotal += w1*w2*scorePair;\r
+ ++uPairCount;\r
+#if TRACE\r
+ Log("%10.2f %6.3f %6.3f %10.2f >%s >%s\n",\r
+ scorePair,\r
+ w1,\r
+ w2,\r
+ scorePair*w1*w2,\r
+ msa1.GetSeqName(uSeqIndex1),\r
+ msa2.GetSeqName(uSeqIndex2));\r
+#endif\r
+ }\r
+ }\r
+ if (0 == uPairCount)\r
+ Quit("0 == uPairCount");\r
+\r
+#if TRACE\r
+ Log("msa1=\n");\r
+ msa1.LogMe();\r
+ Log("msa2=\n");\r
+ msa2.LogMe();\r
+ Log("XP=%g\n", scoreTotal);\r
+#endif\r
+// return scoreTotal / uPairCount;\r
+ return scoreTotal;\r
+ }\r
git://source.jalview.org / jabaws.git / blobdiff