+++ /dev/null
-#include "muscle.h"\r
-#include "msa.h"\r
-#include "profile.h"\r
-#include "objscore.h"\r
-\r
-#if DOUBLE_AFFINE\r
-\r
-#define TRACE 0\r
-#define TEST_SPFAST 0\r
-\r
-static SCORE GapPenalty(unsigned uLength, bool Term, SCORE g, SCORE e)\r
- {\r
- //if (Term)\r
- // {\r
- // switch (g_TermGap)\r
- // {\r
- // case TERMGAP_Full:\r
- // return g + (uLength - 1)*e;\r
-\r
- // case TERMGAP_Half:\r
- // return g/2 + (uLength - 1)*e;\r
-\r
- // case TERMGAP_Ext:\r
- // return uLength*e;\r
- // }\r
- // Quit("Bad termgap");\r
- // }\r
- //else\r
- // return g + (uLength - 1)*e;\r
- //return MINUS_INFINITY;\r
- return g + (uLength - 1)*e;\r
- }\r
-\r
-static SCORE GapPenalty(unsigned uLength, bool Term)\r
- {\r
- SCORE s1 = GapPenalty(uLength, Term, g_scoreGapOpen, g_scoreGapExtend);\r
-#if DOUBLE_AFFINE\r
- SCORE s2 = GapPenalty(uLength, Term, g_scoreGapOpen2, g_scoreGapExtend2);\r
- if (s1 > s2)\r
- return s1;\r
- return s2;\r
-#else\r
- return s1;\r
-#endif\r
- }\r
-\r
-static const MSA *g_ptrMSA1;\r
-static const MSA *g_ptrMSA2;\r
-static unsigned g_uSeqIndex1;\r
-static unsigned g_uSeqIndex2;\r
-\r
-static void LogGap(unsigned uStart, unsigned uEnd, unsigned uGapLength,\r
- bool bNTerm, bool bCTerm)\r
- {\r
- Log("%16.16s ", "");\r
- for (unsigned i = 0; i < uStart; ++i)\r
- Log(" ");\r
- unsigned uMyLength = 0;\r
- for (unsigned i = uStart; i <= uEnd; ++i)\r
- {\r
- bool bGap1 = g_ptrMSA1->IsGap(g_uSeqIndex1, i);\r
- bool bGap2 = g_ptrMSA2->IsGap(g_uSeqIndex2, i);\r
- if (!bGap1 && !bGap2)\r
- Quit("Error -- neither gapping");\r
- if (bGap1 && bGap2)\r
- Log(".");\r
- else\r
- {\r
- ++uMyLength;\r
- Log("-");\r
- }\r
- }\r
- SCORE s = GapPenalty(uGapLength, bNTerm || bCTerm);\r
- Log(" L=%d N%d C%d s=%.3g", uGapLength, bNTerm, bCTerm, s);\r
- Log("\n");\r
- if (uMyLength != uGapLength)\r
- Quit("Lengths differ");\r
-\r
- }\r
-\r
-static SCORE ScoreSeqPair(const MSA &msa1, unsigned uSeqIndex1,\r
- const MSA &msa2, unsigned uSeqIndex2, SCORE *ptrLetters, SCORE *ptrGaps)\r
- {\r
- g_ptrMSA1 = &msa1;\r
- g_ptrMSA2 = &msa2;\r
- g_uSeqIndex1 = uSeqIndex1;\r
- g_uSeqIndex2 = uSeqIndex2;\r
-\r
- const unsigned uColCount = msa1.GetColCount();\r
- const unsigned uColCount2 = msa2.GetColCount();\r
- if (uColCount != uColCount2)\r
- Quit("ScoreSeqPair, different lengths");\r
-\r
-#if TRACE\r
- Log("ScoreSeqPair\n");\r
- Log("%16.16s ", msa1.GetSeqName(uSeqIndex1));\r
- for (unsigned i = 0; i < uColCount; ++i)\r
- Log("%c", msa1.GetChar(uSeqIndex1, i));\r
- Log("\n");\r
- Log("%16.16s ", msa2.GetSeqName(uSeqIndex2));\r
- for (unsigned i = 0; i < uColCount; ++i)\r
- Log("%c", msa1.GetChar(uSeqIndex2, i));\r
- Log("\n");\r
-#endif\r
-\r
- SCORE scoreTotal = 0;\r
-\r
-// Substitution scores\r
- unsigned uFirstLetter1 = uInsane;\r
- unsigned uFirstLetter2 = uInsane;\r
- unsigned uLastLetter1 = uInsane;\r
- unsigned uLastLetter2 = uInsane;\r
- for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)\r
- {\r
- bool bGap1 = msa1.IsGap(uSeqIndex1, uColIndex);\r
- bool bGap2 = msa2.IsGap(uSeqIndex2, uColIndex);\r
- bool bWildcard1 = msa1.IsWildcard(uSeqIndex1, uColIndex);\r
- bool bWildcard2 = msa2.IsWildcard(uSeqIndex2, uColIndex);\r
-\r
- if (!bGap1)\r
- {\r
- if (uInsane == uFirstLetter1)\r
- uFirstLetter1 = uColIndex;\r
- uLastLetter1 = uColIndex;\r
- }\r
- if (!bGap2)\r
- {\r
- if (uInsane == uFirstLetter2)\r
- uFirstLetter2 = uColIndex;\r
- uLastLetter2 = uColIndex;\r
- }\r
-\r
- if (bGap1 || bGap2 || bWildcard1 || bWildcard2)\r
- continue;\r
-\r
- unsigned uLetter1 = msa1.GetLetter(uSeqIndex1, uColIndex);\r
- unsigned uLetter2 = msa2.GetLetter(uSeqIndex2, uColIndex);\r
-\r
- SCORE scoreMatch = (*g_ptrScoreMatrix)[uLetter1][uLetter2];\r
- scoreTotal += scoreMatch;\r
-#if TRACE\r
- Log("%c <-> %c = %7.1f %10.1f\n",\r
- msa1.GetChar(uSeqIndex1, uColIndex),\r
- msa2.GetChar(uSeqIndex2, uColIndex),\r
- scoreMatch,\r
- scoreTotal);\r
-#endif\r
- }\r
- \r
- *ptrLetters = scoreTotal;\r
-\r
-// Gap penalties\r
- unsigned uGapLength = uInsane;\r
- unsigned uGapStartCol = uInsane;\r
- bool bGapping1 = false;\r
- bool bGapping2 = false;\r
-\r
- for (unsigned uColIndex = 0; uColIndex < uColCount; ++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 (bGapping1)\r
- {\r
- if (bGap1)\r
- ++uGapLength;\r
- else\r
- {\r
- bGapping1 = false;\r
- bool bNTerm = (uFirstLetter2 == uGapStartCol);\r
- bool bCTerm = (uLastLetter2 + 1 == uColIndex);\r
- SCORE scoreGap = GapPenalty(uGapLength, bNTerm || bCTerm);\r
- scoreTotal += scoreGap;\r
-#if TRACE\r
- LogGap(uGapStartCol, uColIndex - 1, uGapLength, bNTerm, bCTerm);\r
- Log("GAP %7.1f %10.1f\n",\r
- scoreGap,\r
- scoreTotal);\r
-#endif\r
- }\r
- continue;\r
- }\r
- else\r
- {\r
- if (bGap1)\r
- {\r
- uGapStartCol = uColIndex;\r
- bGapping1 = true;\r
- uGapLength = 1;\r
- continue;\r
- }\r
- }\r
-\r
- if (bGapping2)\r
- {\r
- if (bGap2)\r
- ++uGapLength;\r
- else\r
- {\r
- bGapping2 = false;\r
- bool bNTerm = (uFirstLetter1 == uGapStartCol);\r
- bool bCTerm = (uLastLetter1 + 1 == uColIndex);\r
- SCORE scoreGap = GapPenalty(uGapLength, bNTerm || bCTerm);\r
- scoreTotal += scoreGap;\r
-#if TRACE\r
- LogGap(uGapStartCol, uColIndex - 1, uGapLength, bNTerm, bCTerm);\r
- Log("GAP %7.1f %10.1f\n",\r
- scoreGap,\r
- scoreTotal);\r
-#endif\r
- }\r
- }\r
- else\r
- {\r
- if (bGap2)\r
- {\r
- uGapStartCol = uColIndex;\r
- bGapping2 = true;\r
- uGapLength = 1;\r
- }\r
- }\r
- }\r
-\r
- if (bGapping1 || bGapping2)\r
- {\r
- SCORE scoreGap = GapPenalty(uGapLength, true);\r
- scoreTotal += scoreGap;\r
-#if TRACE\r
- LogGap(uGapStartCol, uColCount - 1, uGapLength, false, true);\r
- Log("GAP %7.1f %10.1f\n",\r
- scoreGap,\r
- scoreTotal);\r
-#endif\r
- }\r
- *ptrGaps = scoreTotal - *ptrLetters;\r
- return scoreTotal;\r
- }\r
-\r
-// The usual sum-of-pairs objective score: sum the score\r
-// of the alignment of each pair of sequences.\r
-SCORE ObjScoreDA(const MSA &msa, SCORE *ptrLetters, SCORE *ptrGaps)\r
- {\r
- const unsigned uSeqCount = msa.GetSeqCount();\r
- SCORE scoreTotal = 0;\r
- unsigned uPairCount = 0;\r
-#if TRACE\r
- msa.LogMe();\r
- Log(" Score Weight Weight Total\n");\r
- Log("---------- ------ ------ ----------\n");\r
-#endif\r
- SCORE TotalLetters = 0;\r
- SCORE TotalGaps = 0;\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
- SCORE Letters;\r
- SCORE Gaps;\r
- SCORE scorePair = ScoreSeqPair(msa, uSeqIndex1, msa, uSeqIndex2,\r
- &Letters, &Gaps);\r
- scoreTotal += w1*w2*scorePair;\r
- TotalLetters += w1*w2*Letters;\r
- TotalGaps += w1*w2*Gaps;\r
- ++uPairCount;\r
-#if TRACE\r
- Log("%10.2f %6.3f %6.3f %10.2f %d=%s %d=%s\n",\r
- scorePair,\r
- w1,\r
- w2,\r
- scorePair*w1*w2,\r
- uSeqIndex1,\r
- msa.GetSeqName(uSeqIndex1),\r
- uSeqIndex2,\r
- msa.GetSeqName(uSeqIndex2));\r
-#endif\r
- }\r
- }\r
- *ptrLetters = TotalLetters;\r
- *ptrGaps = TotalGaps;\r
- return scoreTotal;\r
- }\r
-\r
-#endif // DOUBLE_AFFINE\r
git://source.jalview.org / jabaws.git / blobdiff