website/archive/binaries/mac/src/muscle/diffobjscore.cpp

   1 #include "muscle.h"\r
   2 #include "msa.h"\r
   3 #include "objscore.h"\r
   4 #include "profile.h"\r
   5 \r
   6 #define TRACE                           0\r
   7 #define COMPARE_3_52            0\r
   8 #define BRUTE_LETTERS           0\r
   9 \r
  10 static SCORE ScoreColLetters(const MSA &msa, unsigned uColIndex)\r
  11         {\r
  12         SCOREMATRIX &Mx = *g_ptrScoreMatrix;\r
  13         const unsigned uSeqCount = msa.GetSeqCount();\r
  14 \r
  15 #if     BRUTE_LETTERS\r
  16         SCORE BruteScore = 0;\r
  17         for (unsigned uSeqIndex1 = 0; uSeqIndex1 < uSeqCount; ++uSeqIndex1)\r
  18                 {\r
  19                 unsigned uLetter1 = msa.GetLetterEx(uSeqIndex1, uColIndex);\r
  20                 if (uLetter1 >= g_AlphaSize)\r
  21                         continue;\r
  22                 WEIGHT w1 = msa.GetSeqWeight(uSeqIndex1);\r
  23                 for (unsigned uSeqIndex2 = uSeqIndex1 + 1; uSeqIndex2 < uSeqCount; ++uSeqIndex2)\r
  24                         {\r
  25                         unsigned uLetter2 = msa.GetLetterEx(uSeqIndex2, uColIndex);\r
  26                         if (uLetter2 >= g_AlphaSize)\r
  27                                 continue;\r
  28                         WEIGHT w2 = msa.GetSeqWeight(uSeqIndex2);\r
  29                         BruteScore += w1*w2*Mx[uLetter1][uLetter2];\r
  30                         }\r
  31                 }\r
  32 #endif\r
  33         \r
  34         double N = 0;\r
  35         for (unsigned uSeqIndex1 = 0; uSeqIndex1 < uSeqCount; ++uSeqIndex1)\r
  36                 {\r
  37                 WEIGHT w = msa.GetSeqWeight(uSeqIndex1);\r
  38                 N += w;\r
  39                 }\r
  40         if (N <= 0)\r
  41                 return 0;\r
  42 \r
  43         FCOUNT Freqs[20];\r
  44         memset(Freqs, 0, sizeof(Freqs));\r
  45         SCORE Score = 0;\r
  46         for (unsigned uSeqIndex1 = 0; uSeqIndex1 < uSeqCount; ++uSeqIndex1)\r
  47                 {\r
  48                 unsigned uLetter = msa.GetLetterEx(uSeqIndex1, uColIndex);\r
  49                 if (uLetter >= g_AlphaSize)\r
  50                         continue;\r
  51                 WEIGHT w = msa.GetSeqWeight(uSeqIndex1);\r
  52                 Freqs[uLetter] += w;\r
  53                 Score -= w*w*Mx[uLetter][uLetter];\r
  54                 }\r
  55 \r
  56         for (unsigned uLetter1 = 0; uLetter1 < g_AlphaSize; ++uLetter1)\r
  57                 {\r
  58                 const FCOUNT f1 = Freqs[uLetter1];\r
  59                 Score += f1*f1*Mx[uLetter1][uLetter1];\r
  60                 for (unsigned uLetter2 = uLetter1 + 1; uLetter2 < g_AlphaSize; ++uLetter2)\r
  61                         {\r
  62                         const FCOUNT f2 = Freqs[uLetter2];\r
  63                         Score += 2*f1*f2*Mx[uLetter1][uLetter2];\r
  64                         }\r
  65                 }\r
  66         Score /= 2;\r
  67 #if     BRUTE_LETTERS\r
  68         assert(BTEq(BruteScore, Score));\r
  69 #endif\r
  70         return Score;\r
  71         }\r
  72 \r
  73 static SCORE ScoreLetters(const MSA &msa, const unsigned Edges[],\r
  74   unsigned uEdgeCount)\r
  75         {\r
  76         const unsigned uSeqCount = msa.GetSeqCount();\r
  77         const unsigned uColCount = msa.GetColCount();\r
  78 \r
  79 // Letters\r
  80         SCORE Score = 0;\r
  81         for (unsigned uEdgeIndex = 0; uEdgeIndex < uEdgeCount; ++uEdgeIndex)\r
  82                 {\r
  83                 const unsigned uColIndex = Edges[uEdgeIndex];\r
  84                 assert(uColIndex < uColCount);\r
  85                 Score += ScoreColLetters(msa, uColIndex);\r
  86                 }\r
  87         return Score;\r
  88         }\r
  89 \r
  90 void GetLetterScores(const MSA &msa, SCORE Scores[])\r
  91         {\r
  92         const unsigned uColCount = msa.GetColCount();\r
  93         for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)\r
  94                 Scores[uColIndex] = ScoreColLetters(msa, uColIndex);\r
  95         }\r
  96 \r
  97 SCORE DiffObjScore(\r
  98   const MSA &msa1, const PWPath &Path1, const unsigned Edges1[], unsigned uEdgeCount1, \r
  99   const MSA &msa2, const PWPath &Path2, const unsigned Edges2[], unsigned uEdgeCount2)\r
 100         {\r
 101 #if     TRACE\r
 102         {\r
 103         Log("============DiffObjScore===========\n");\r
 104         Log("msa1:\n");\r
 105         msa1.LogMe();\r
 106         Log("\n");\r
 107         Log("Cols1: ");\r
 108         for (unsigned i = 0; i < uEdgeCount1; ++i)\r
 109                 Log(" %u", Edges1[i]);\r
 110         Log("\n\n");\r
 111         Log("msa2:\n");\r
 112         msa2.LogMe();\r
 113         Log("Cols2: ");\r
 114         for (unsigned i = 0; i < uEdgeCount2; ++i)\r
 115                 Log(" %u", Edges2[i]);\r
 116         Log("\n\n");\r
 117         }\r
 118 #endif\r
 119 \r
 120 #if     COMPARE_3_52\r
 121         extern SCORE g_SPScoreLetters;\r
 122         extern SCORE g_SPScoreGaps;\r
 123         SCORE SP1 = ObjScoreSP(msa1);\r
 124         SCORE SPLetters1 = g_SPScoreLetters;\r
 125         SCORE SPGaps1 = g_SPScoreGaps;\r
 126 \r
 127         SCORE SP2 = ObjScoreSP(msa2);\r
 128         SCORE SPLetters2 = g_SPScoreLetters;\r
 129         SCORE SPGaps2 = g_SPScoreGaps;\r
 130         SCORE SPDiffLetters = SPLetters2 - SPLetters1;\r
 131         SCORE SPDiffGaps = SPGaps2 - SPGaps1;\r
 132         SCORE SPDiff = SPDiffLetters + SPDiffGaps;\r
 133 #endif\r
 134 \r
 135         SCORE Letters1 = ScoreLetters(msa1, Edges1, uEdgeCount1);\r
 136         SCORE Letters2 = ScoreLetters(msa2, Edges2, uEdgeCount2);\r
 137 \r
 138         SCORE Gaps1 = ScoreGaps(msa1, Edges1, uEdgeCount1);\r
 139         SCORE Gaps2 = ScoreGaps(msa2, Edges2, uEdgeCount2);\r
 140 \r
 141         SCORE DiffLetters = Letters2 - Letters1;\r
 142         SCORE DiffGaps = Gaps2 - Gaps1;\r
 143         SCORE Diff = DiffLetters + DiffGaps;\r
 144 \r
 145 #if     COMPARE_3_52\r
 146         Log("ObjScoreSP    Letters1=%.4g  Letters2=%.4g  DiffLetters=%.4g\n",\r
 147           SPLetters1, SPLetters2, SPDiffLetters);\r
 148 \r
 149         Log("DiffObjScore  Letters1=%.4g  Letters2=%.4g  DiffLetters=%.4g\n",\r
 150           Letters1, Letters2, DiffLetters);\r
 151 \r
 152         Log("ObjScoreSP    Gaps1=%.4g  Gaps2=%.4g  DiffGaps=%.4g\n",\r
 153           SPGaps1, SPGaps2, SPDiffGaps);\r
 154 \r
 155         Log("DiffObjScore  Gaps1=%.4g  Gaps2=%.4g  DiffGaps=%.4g\n",\r
 156           Gaps1, Gaps2, DiffGaps);\r
 157 \r
 158         Log("SP diff=%.4g DiffObjScore Diff=%.4g\n", SPDiff, Diff);\r
 159 #endif\r
 160 \r
 161         return Diff;\r
 162         }\r