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

   1 #include "muscle.h"\r
   2 #include "msa.h"\r
   3 #include "profile.h"\r
   4 #include "objscore.h"\r
   5 \r
   6 #if     DOUBLE_AFFINE\r
   7 \r
   8 #define TRACE                   0\r
   9 #define TEST_SPFAST             0\r
  10 \r
  11 static SCORE GapPenalty(unsigned uLength, bool Term, SCORE g, SCORE e)\r
  12         {\r
  13         //if (Term)\r
  14         //      {\r
  15         //      switch (g_TermGap)\r
  16         //              {\r
  17         //      case TERMGAP_Full:\r
  18         //              return g + (uLength - 1)*e;\r
  19 \r
  20         //      case TERMGAP_Half:\r
  21         //              return g/2 + (uLength - 1)*e;\r
  22 \r
  23         //      case TERMGAP_Ext:\r
  24         //              return uLength*e;\r
  25         //              }\r
  26         //      Quit("Bad termgap");\r
  27         //      }\r
  28         //else\r
  29         //      return g + (uLength - 1)*e;\r
  30         //return MINUS_INFINITY;\r
  31         return g + (uLength - 1)*e;\r
  32         }\r
  33 \r
  34 static SCORE GapPenalty(unsigned uLength, bool Term)\r
  35         {\r
  36         SCORE s1 = GapPenalty(uLength, Term, g_scoreGapOpen, g_scoreGapExtend);\r
  37 #if     DOUBLE_AFFINE\r
  38         SCORE s2 = GapPenalty(uLength, Term, g_scoreGapOpen2, g_scoreGapExtend2);\r
  39         if (s1 > s2)\r
  40                 return s1;\r
  41         return s2;\r
  42 #else\r
  43         return s1;\r
  44 #endif\r
  45         }\r
  46 \r
  47 static const MSA *g_ptrMSA1;\r
  48 static const MSA *g_ptrMSA2;\r
  49 static unsigned g_uSeqIndex1;\r
  50 static unsigned g_uSeqIndex2;\r
  51 \r
  52 static void LogGap(unsigned uStart, unsigned uEnd, unsigned uGapLength,\r
  53   bool bNTerm, bool bCTerm)\r
  54         {\r
  55         Log("%16.16s  ", "");\r
  56         for (unsigned i = 0; i < uStart; ++i)\r
  57                 Log(" ");\r
  58         unsigned uMyLength = 0;\r
  59         for (unsigned i = uStart; i <= uEnd; ++i)\r
  60                 {\r
  61                 bool bGap1 = g_ptrMSA1->IsGap(g_uSeqIndex1, i);\r
  62                 bool bGap2 = g_ptrMSA2->IsGap(g_uSeqIndex2, i);\r
  63                 if (!bGap1 && !bGap2)\r
  64                         Quit("Error -- neither gapping");\r
  65                 if (bGap1 && bGap2)\r
  66                         Log(".");\r
  67                 else\r
  68                         {\r
  69                         ++uMyLength;\r
  70                         Log("-");\r
  71                         }\r
  72                 }\r
  73         SCORE s = GapPenalty(uGapLength, bNTerm || bCTerm);\r
  74         Log(" L=%d N%d C%d s=%.3g", uGapLength, bNTerm, bCTerm, s);\r
  75         Log("\n");\r
  76         if (uMyLength != uGapLength)\r
  77                 Quit("Lengths differ");\r
  78 \r
  79         }\r
  80 \r
  81 static SCORE ScoreSeqPair(const MSA &msa1, unsigned uSeqIndex1,\r
  82   const MSA &msa2, unsigned uSeqIndex2, SCORE *ptrLetters, SCORE *ptrGaps)\r
  83         {\r
  84         g_ptrMSA1 = &msa1;\r
  85         g_ptrMSA2 = &msa2;\r
  86         g_uSeqIndex1 = uSeqIndex1;\r
  87         g_uSeqIndex2 = uSeqIndex2;\r
  88 \r
  89         const unsigned uColCount = msa1.GetColCount();\r
  90         const unsigned uColCount2 = msa2.GetColCount();\r
  91         if (uColCount != uColCount2)\r
  92                 Quit("ScoreSeqPair, different lengths");\r
  93 \r
  94 #if     TRACE\r
  95         Log("ScoreSeqPair\n");\r
  96         Log("%16.16s  ", msa1.GetSeqName(uSeqIndex1));\r
  97         for (unsigned i = 0; i < uColCount; ++i)\r
  98                 Log("%c", msa1.GetChar(uSeqIndex1, i));\r
  99         Log("\n");\r
 100         Log("%16.16s  ", msa2.GetSeqName(uSeqIndex2));\r
 101         for (unsigned i = 0; i < uColCount; ++i)\r
 102                 Log("%c", msa1.GetChar(uSeqIndex2, i));\r
 103         Log("\n");\r
 104 #endif\r
 105 \r
 106         SCORE scoreTotal = 0;\r
 107 \r
 108 // Substitution scores\r
 109         unsigned uFirstLetter1 = uInsane;\r
 110         unsigned uFirstLetter2 = uInsane;\r
 111         unsigned uLastLetter1 = uInsane;\r
 112         unsigned uLastLetter2 = uInsane;\r
 113         for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)\r
 114                 {\r
 115                 bool bGap1 = msa1.IsGap(uSeqIndex1, uColIndex);\r
 116                 bool bGap2 = msa2.IsGap(uSeqIndex2, uColIndex);\r
 117                 bool bWildcard1 = msa1.IsWildcard(uSeqIndex1, uColIndex);\r
 118                 bool bWildcard2 = msa2.IsWildcard(uSeqIndex2, uColIndex);\r
 119 \r
 120                 if (!bGap1)\r
 121                         {\r
 122                         if (uInsane == uFirstLetter1)\r
 123                                 uFirstLetter1 = uColIndex;\r
 124                         uLastLetter1 = uColIndex;\r
 125                         }\r
 126                 if (!bGap2)\r
 127                         {\r
 128                         if (uInsane == uFirstLetter2)\r
 129                                 uFirstLetter2 = uColIndex;\r
 130                         uLastLetter2 = uColIndex;\r
 131                         }\r
 132 \r
 133                 if (bGap1 || bGap2 || bWildcard1 || bWildcard2)\r
 134                         continue;\r
 135 \r
 136                 unsigned uLetter1 = msa1.GetLetter(uSeqIndex1, uColIndex);\r
 137                 unsigned uLetter2 = msa2.GetLetter(uSeqIndex2, uColIndex);\r
 138 \r
 139                 SCORE scoreMatch = (*g_ptrScoreMatrix)[uLetter1][uLetter2];\r
 140                 scoreTotal += scoreMatch;\r
 141 #if     TRACE\r
 142                 Log("%c <-> %c = %7.1f  %10.1f\n",\r
 143                   msa1.GetChar(uSeqIndex1, uColIndex),\r
 144                   msa2.GetChar(uSeqIndex2, uColIndex),\r
 145                   scoreMatch,\r
 146                   scoreTotal);\r
 147 #endif\r
 148                 }\r
 149         \r
 150         *ptrLetters = scoreTotal;\r
 151 \r
 152 // Gap penalties\r
 153         unsigned uGapLength = uInsane;\r
 154         unsigned uGapStartCol = uInsane;\r
 155         bool bGapping1 = false;\r
 156         bool bGapping2 = false;\r
 157 \r
 158         for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)\r
 159                 {\r
 160                 bool bGap1 = msa1.IsGap(uSeqIndex1, uColIndex);\r
 161                 bool bGap2 = msa2.IsGap(uSeqIndex2, uColIndex);\r
 162 \r
 163                 if (bGap1 && bGap2)\r
 164                         continue;\r
 165 \r
 166                 if (bGapping1)\r
 167                         {\r
 168                         if (bGap1)\r
 169                                 ++uGapLength;\r
 170                         else\r
 171                                 {\r
 172                                 bGapping1 = false;\r
 173                                 bool bNTerm = (uFirstLetter2 == uGapStartCol);\r
 174                                 bool bCTerm = (uLastLetter2 + 1 == uColIndex);\r
 175                                 SCORE scoreGap = GapPenalty(uGapLength, bNTerm || bCTerm);\r
 176                                 scoreTotal += scoreGap;\r
 177 #if     TRACE\r
 178                                 LogGap(uGapStartCol, uColIndex - 1, uGapLength, bNTerm, bCTerm);\r
 179                                 Log("GAP         %7.1f  %10.1f\n",\r
 180                                   scoreGap,\r
 181                                   scoreTotal);\r
 182 #endif\r
 183                                 }\r
 184                         continue;\r
 185                         }\r
 186                 else\r
 187                         {\r
 188                         if (bGap1)\r
 189                                 {\r
 190                                 uGapStartCol = uColIndex;\r
 191                                 bGapping1 = true;\r
 192                                 uGapLength = 1;\r
 193                                 continue;\r
 194                                 }\r
 195                         }\r
 196 \r
 197                 if (bGapping2)\r
 198                         {\r
 199                         if (bGap2)\r
 200                                 ++uGapLength;\r
 201                         else\r
 202                                 {\r
 203                                 bGapping2 = false;\r
 204                                 bool bNTerm = (uFirstLetter1 == uGapStartCol);\r
 205                                 bool bCTerm = (uLastLetter1 + 1 == uColIndex);\r
 206                                 SCORE scoreGap = GapPenalty(uGapLength, bNTerm || bCTerm);\r
 207                                 scoreTotal += scoreGap;\r
 208 #if     TRACE\r
 209                                 LogGap(uGapStartCol, uColIndex - 1, uGapLength, bNTerm, bCTerm);\r
 210                                 Log("GAP         %7.1f  %10.1f\n",\r
 211                                   scoreGap,\r
 212                                   scoreTotal);\r
 213 #endif\r
 214                                 }\r
 215                         }\r
 216                 else\r
 217                         {\r
 218                         if (bGap2)\r
 219                                 {\r
 220                                 uGapStartCol = uColIndex;\r
 221                                 bGapping2 = true;\r
 222                                 uGapLength = 1;\r
 223                                 }\r
 224                         }\r
 225                 }\r
 226 \r
 227         if (bGapping1 || bGapping2)\r
 228                 {\r
 229                 SCORE scoreGap = GapPenalty(uGapLength, true);\r
 230                 scoreTotal += scoreGap;\r
 231 #if     TRACE\r
 232                 LogGap(uGapStartCol, uColCount - 1, uGapLength, false, true);\r
 233                 Log("GAP         %7.1f  %10.1f\n",\r
 234                   scoreGap,\r
 235                   scoreTotal);\r
 236 #endif\r
 237                 }\r
 238         *ptrGaps = scoreTotal - *ptrLetters;\r
 239         return scoreTotal;\r
 240         }\r
 241 \r
 242 // The usual sum-of-pairs objective score: sum the score\r
 243 // of the alignment of each pair of sequences.\r
 244 SCORE ObjScoreDA(const MSA &msa, SCORE *ptrLetters, SCORE *ptrGaps)\r
 245         {\r
 246         const unsigned uSeqCount = msa.GetSeqCount();\r
 247         SCORE scoreTotal = 0;\r
 248         unsigned uPairCount = 0;\r
 249 #if     TRACE\r
 250         msa.LogMe();\r
 251         Log("     Score  Weight  Weight       Total\n");\r
 252         Log("----------  ------  ------  ----------\n");\r
 253 #endif\r
 254         SCORE TotalLetters = 0;\r
 255         SCORE TotalGaps = 0;\r
 256         for (unsigned uSeqIndex1 = 0; uSeqIndex1 < uSeqCount; ++uSeqIndex1)\r
 257                 {\r
 258                 const WEIGHT w1 = msa.GetSeqWeight(uSeqIndex1);\r
 259                 for (unsigned uSeqIndex2 = uSeqIndex1 + 1; uSeqIndex2 < uSeqCount; ++uSeqIndex2)\r
 260                         {\r
 261                         const WEIGHT w2 = msa.GetSeqWeight(uSeqIndex2);\r
 262                         const WEIGHT w = w1*w2;\r
 263                         SCORE Letters;\r
 264                         SCORE Gaps;\r
 265                         SCORE scorePair = ScoreSeqPair(msa, uSeqIndex1, msa, uSeqIndex2,\r
 266                           &Letters, &Gaps);\r
 267                         scoreTotal += w1*w2*scorePair;\r
 268                         TotalLetters += w1*w2*Letters;\r
 269                         TotalGaps += w1*w2*Gaps;\r
 270                         ++uPairCount;\r
 271 #if     TRACE\r
 272                         Log("%10.2f  %6.3f  %6.3f  %10.2f  %d=%s %d=%s\n",\r
 273                           scorePair,\r
 274                           w1,\r
 275                           w2,\r
 276                           scorePair*w1*w2,\r
 277                           uSeqIndex1,\r
 278                           msa.GetSeqName(uSeqIndex1),\r
 279                           uSeqIndex2,\r
 280                           msa.GetSeqName(uSeqIndex2));\r
 281 #endif\r
 282                         }\r
 283                 }\r
 284         *ptrLetters = TotalLetters;\r
 285         *ptrGaps = TotalGaps;\r
 286         return scoreTotal;\r
 287         }\r
 288 \r
 289 #endif  // DOUBLE_AFFINE\r