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

   1 #include "muscle.h"\r
   2 #include <math.h>\r
   3 #include "pwpath.h"\r
   4 #include "profile.h"\r
   5 #include <stdio.h>\r
   6 \r
   7 #define TRACE   0\r
   8 \r
   9 #if     1 // SINGLE_AFFINE\r
  10 \r
  11 extern bool g_bKeepSimpleDP;\r
  12 extern SCORE *g_DPM;\r
  13 extern SCORE *g_DPD;\r
  14 extern SCORE *g_DPI;\r
  15 extern char *g_TBM;\r
  16 extern char *g_TBD;\r
  17 extern char *g_TBI;\r
  18 \r
  19 static const char *LocalScoreToStr(SCORE s)\r
  20         {\r
  21         static char str[16];\r
  22         if (s < -100000)\r
  23                 return "     *";\r
  24         sprintf(str, "%6.1f", s);\r
  25         return str;\r
  26         }\r
  27 \r
  28 static void ListTB(const char *TBM_, const ProfPos *PA, const ProfPos *PB,\r
  29   unsigned uPrefixCountA, unsigned uPrefixCountB)\r
  30         {\r
  31         Log("        ");\r
  32         for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
  33                 {\r
  34                 char c = ' ';\r
  35                 if (uPrefixLengthB > 0)\r
  36                         c = ConsensusChar(PB[uPrefixLengthB - 1]);\r
  37                 Log(" %4u:%c", uPrefixLengthB, c);\r
  38                 }\r
  39         Log("\n");\r
  40         for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)\r
  41                 {\r
  42                 char c = ' ';\r
  43                 if (uPrefixLengthA > 0)\r
  44                         c = ConsensusChar(PA[uPrefixLengthA - 1]);\r
  45                 Log("%4u:%c  ", uPrefixLengthA, c);\r
  46                 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
  47                         Log(" %6c", TBM(uPrefixLengthA, uPrefixLengthB));\r
  48                 Log("\n");\r
  49                 }\r
  50         }\r
  51 \r
  52 static void ListDP(const SCORE *DPM_, const ProfPos *PA, const ProfPos *PB,\r
  53   unsigned uPrefixCountA, unsigned uPrefixCountB)\r
  54         {\r
  55         Log("        ");\r
  56         for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
  57                 {\r
  58                 char c = ' ';\r
  59                 if (uPrefixLengthB > 0)\r
  60                         c = ConsensusChar(PB[uPrefixLengthB - 1]);\r
  61                 Log(" %4u:%c", uPrefixLengthB, c);\r
  62                 }\r
  63         Log("\n");\r
  64         for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)\r
  65                 {\r
  66                 char c = ' ';\r
  67                 if (uPrefixLengthA > 0)\r
  68                         c = ConsensusChar(PA[uPrefixLengthA - 1]);\r
  69                 Log("%4u:%c  ", uPrefixLengthA, c);\r
  70                 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
  71                         Log(" %s", LocalScoreToStr(DPM(uPrefixLengthA, uPrefixLengthB)));\r
  72                 Log("\n");\r
  73                 }\r
  74         }\r
  75 \r
  76 SCORE GlobalAlignSimple(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,\r
  77   unsigned uLengthB, PWPath &Path)\r
  78         {\r
  79         assert(uLengthB > 0 && uLengthA > 0);\r
  80 \r
  81         SetTermGaps(PA, uLengthA);\r
  82         SetTermGaps(PB, uLengthB);\r
  83 \r
  84         const unsigned uPrefixCountA = uLengthA + 1;\r
  85         const unsigned uPrefixCountB = uLengthB + 1;\r
  86 \r
  87 // Allocate DP matrices\r
  88         const size_t LM = uPrefixCountA*uPrefixCountB;\r
  89         SCORE *DPL_ = new SCORE[LM];\r
  90         SCORE *DPM_ = new SCORE[LM];\r
  91         SCORE *DPD_ = new SCORE[LM];\r
  92         SCORE *DPI_ = new SCORE[LM];\r
  93 \r
  94         char *TBM_ = new char[LM];\r
  95         char *TBD_ = new char[LM];\r
  96         char *TBI_ = new char[LM];\r
  97 \r
  98         memset(TBM_, '?', LM);\r
  99         memset(TBD_, '?', LM);\r
 100         memset(TBI_, '?', LM);\r
 101 \r
 102         DPM(0, 0) = 0;\r
 103         DPD(0, 0) = MINUS_INFINITY;\r
 104         DPI(0, 0) = MINUS_INFINITY;\r
 105 \r
 106         DPM(1, 0) = MINUS_INFINITY;\r
 107         DPD(1, 0) = PA[0].m_scoreGapOpen;\r
 108         TBD(1, 0) = 'D';\r
 109         DPI(1, 0) = MINUS_INFINITY;\r
 110 \r
 111         DPM(0, 1) = MINUS_INFINITY;\r
 112         DPD(0, 1) = MINUS_INFINITY;\r
 113         DPI(0, 1) = PB[0].m_scoreGapOpen;\r
 114         TBI(0, 1) = 'I';\r
 115 \r
 116 // Empty prefix of B is special case\r
 117         for (unsigned uPrefixLengthA = 2; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)\r
 118                 {\r
 119         // M=LetterA+LetterB, impossible with empty prefix\r
 120                 DPM(uPrefixLengthA, 0) = MINUS_INFINITY;\r
 121 \r
 122         // D=LetterA+GapB\r
 123                 DPD(uPrefixLengthA, 0) = DPD(uPrefixLengthA - 1, 0) + g_scoreGapExtend;\r
 124                 TBD(uPrefixLengthA, 0) = 'D';\r
 125 \r
 126         // I=GapA+LetterB, impossible with empty prefix\r
 127                 DPI(uPrefixLengthA, 0) = MINUS_INFINITY;\r
 128                 }\r
 129 \r
 130 // Empty prefix of A is special case\r
 131         for (unsigned uPrefixLengthB = 2; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
 132                 {\r
 133         // M=LetterA+LetterB, impossible with empty prefix\r
 134                 DPM(0, uPrefixLengthB) = MINUS_INFINITY;\r
 135 \r
 136         // D=LetterA+GapB, impossible with empty prefix\r
 137                 DPD(0, uPrefixLengthB) = MINUS_INFINITY;\r
 138 \r
 139         // I=GapA+LetterB\r
 140                 DPI(0, uPrefixLengthB) = DPI(0, uPrefixLengthB - 1) + g_scoreGapExtend;\r
 141                 TBI(0, uPrefixLengthB) = 'I';\r
 142                 }\r
 143 \r
 144 // Special case to agree with NWFast, no D-I transitions so...\r
 145         DPD(uLengthA, 0) = MINUS_INFINITY;\r
 146 //      DPI(0, uLengthB) = MINUS_INFINITY;\r
 147 \r
 148 // ============\r
 149 // Main DP loop\r
 150 // ============\r
 151         SCORE scoreGapCloseB = MINUS_INFINITY;\r
 152         for (unsigned uPrefixLengthB = 1; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
 153                 {\r
 154                 const ProfPos &PPB = PB[uPrefixLengthB - 1];\r
 155 \r
 156                 SCORE scoreGapCloseA = MINUS_INFINITY;\r
 157                 for (unsigned uPrefixLengthA = 1; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)\r
 158                         {\r
 159                         const ProfPos &PPA = PA[uPrefixLengthA - 1];\r
 160 \r
 161                         {\r
 162                 // Match M=LetterA+LetterB\r
 163                         SCORE scoreLL = ScoreProfPos2(PPA, PPB);\r
 164                         DPL(uPrefixLengthA, uPrefixLengthB) = scoreLL;\r
 165 \r
 166                         SCORE scoreMM = DPM(uPrefixLengthA-1, uPrefixLengthB-1);\r
 167                         SCORE scoreDM = DPD(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseA;\r
 168                         SCORE scoreIM = DPI(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseB;\r
 169 \r
 170                         SCORE scoreBest;\r
 171                         if (scoreMM >= scoreDM && scoreMM >= scoreIM)\r
 172                                 {\r
 173                                 scoreBest = scoreMM;\r
 174                                 TBM(uPrefixLengthA, uPrefixLengthB) = 'M';\r
 175                                 }\r
 176                         else if (scoreDM >= scoreMM && scoreDM >= scoreIM)\r
 177                                 {\r
 178                                 scoreBest = scoreDM;\r
 179                                 TBM(uPrefixLengthA, uPrefixLengthB) = 'D';\r
 180                                 }\r
 181                         else \r
 182                                 {\r
 183                                 assert(scoreIM >= scoreMM && scoreIM >= scoreDM);\r
 184                                 scoreBest = scoreIM;\r
 185                                 TBM(uPrefixLengthA, uPrefixLengthB) = 'I';\r
 186                                 }\r
 187                         DPM(uPrefixLengthA, uPrefixLengthB) = scoreBest + scoreLL;\r
 188                         }\r
 189 \r
 190                         {\r
 191                 // Delete D=LetterA+GapB\r
 192                         SCORE scoreMD = DPM(uPrefixLengthA-1, uPrefixLengthB) +\r
 193                           PA[uPrefixLengthA-1].m_scoreGapOpen;\r
 194                         SCORE scoreDD = DPD(uPrefixLengthA-1, uPrefixLengthB) + g_scoreGapExtend;\r
 195 \r
 196                         SCORE scoreBest;\r
 197                         if (scoreMD >= scoreDD)\r
 198                                 {\r
 199                                 scoreBest = scoreMD;\r
 200                                 TBD(uPrefixLengthA, uPrefixLengthB) = 'M';\r
 201                                 }\r
 202                         else\r
 203                                 {\r
 204                                 assert(scoreDD >= scoreMD);\r
 205                                 scoreBest = scoreDD;\r
 206                                 TBD(uPrefixLengthA, uPrefixLengthB) = 'D';\r
 207                                 }\r
 208                         DPD(uPrefixLengthA, uPrefixLengthB) = scoreBest;\r
 209                         }\r
 210 \r
 211                 // Insert I=GapA+LetterB\r
 212                         {\r
 213                         SCORE scoreMI = DPM(uPrefixLengthA, uPrefixLengthB-1) +\r
 214                           PB[uPrefixLengthB - 1].m_scoreGapOpen;\r
 215                         SCORE scoreII = DPI(uPrefixLengthA, uPrefixLengthB-1) + g_scoreGapExtend;\r
 216 \r
 217                         SCORE scoreBest;\r
 218                         if (scoreMI >= scoreII)\r
 219                                 {\r
 220                                 scoreBest = scoreMI;\r
 221                                 TBI(uPrefixLengthA, uPrefixLengthB) = 'M';\r
 222                                 }\r
 223                         else \r
 224                                 {\r
 225                                 assert(scoreII > scoreMI);\r
 226                                 scoreBest = scoreII;\r
 227                                 TBI(uPrefixLengthA, uPrefixLengthB) = 'I';\r
 228                                 }\r
 229                         DPI(uPrefixLengthA, uPrefixLengthB) = scoreBest;\r
 230                         }\r
 231 \r
 232                         scoreGapCloseA = PPA.m_scoreGapClose;\r
 233                         }\r
 234                 scoreGapCloseB = PPB.m_scoreGapClose;\r
 235                 }\r
 236 \r
 237 #if TRACE\r
 238         Log("\n");\r
 239         Log("Simple DPL:\n");\r
 240         ListDP(DPL_, PA, PB, uPrefixCountA, uPrefixCountB);\r
 241         Log("\n");\r
 242         Log("Simple DPM:\n");\r
 243         ListDP(DPM_, PA, PB, uPrefixCountA, uPrefixCountB);\r
 244         Log("\n");\r
 245         Log("Simple DPD:\n");\r
 246         ListDP(DPD_, PA, PB, uPrefixCountA, uPrefixCountB);\r
 247         Log("\n");\r
 248         Log("Simple DPI:\n");\r
 249         ListDP(DPI_, PA, PB, uPrefixCountA, uPrefixCountB);\r
 250         Log("\n");\r
 251         Log("Simple TBM:\n");\r
 252         ListTB(TBM_, PA, PB, uPrefixCountA, uPrefixCountB);\r
 253         Log("\n");\r
 254         Log("Simple TBD:\n");\r
 255         ListTB(TBD_, PA, PB, uPrefixCountA, uPrefixCountB);\r
 256         Log("\n");\r
 257         Log("Simple TBI:\n");\r
 258         ListTB(TBI_, PA, PB, uPrefixCountA, uPrefixCountB);\r
 259 #endif\r
 260 \r
 261 // Trace-back\r
 262 // ==========\r
 263         Path.Clear();\r
 264 \r
 265 // Find last edge\r
 266         SCORE M = DPM(uLengthA, uLengthB);\r
 267         SCORE D = DPD(uLengthA, uLengthB) + PA[uLengthA-1].m_scoreGapClose;\r
 268         SCORE I = DPI(uLengthA, uLengthB) + PB[uLengthB-1].m_scoreGapClose;\r
 269         char cEdgeType = '?';\r
 270 \r
 271         SCORE BestScore = MINUS_INFINITY;\r
 272         if (M >= D && M >= I)\r
 273                 {\r
 274                 cEdgeType = 'M';\r
 275                 BestScore = M;\r
 276                 }\r
 277         else if (D >= M && D >= I)\r
 278                 {\r
 279                 cEdgeType = 'D';\r
 280                 BestScore = D;\r
 281                 }\r
 282         else \r
 283                 {\r
 284                 assert(I >= M && I >= D);\r
 285                 cEdgeType = 'I';\r
 286                 BestScore = I;\r
 287                 }\r
 288 \r
 289 #if     TRACE\r
 290         Log("Simple: MAB=%.4g DAB=%.4g IAB=%.4g best=%c\n", M, D, I, cEdgeType);\r
 291 #endif\r
 292 \r
 293         unsigned PLA = uLengthA;\r
 294         unsigned PLB = uLengthB;\r
 295         for (;;)\r
 296                 {\r
 297                 PWEdge Edge;\r
 298                 Edge.cType = cEdgeType;\r
 299                 Edge.uPrefixLengthA = PLA;\r
 300                 Edge.uPrefixLengthB = PLB;\r
 301 #if     TRACE\r
 302                 Log("Prepend %c%d.%d\n", Edge.cType, PLA, PLB);\r
 303 #endif\r
 304                 Path.PrependEdge(Edge);\r
 305 \r
 306                 switch (cEdgeType)\r
 307                         {\r
 308                 case 'M':\r
 309                         assert(PLA > 0);\r
 310                         assert(PLB > 0);\r
 311                         cEdgeType = TBM(PLA, PLB);\r
 312                         --PLA;\r
 313                         --PLB;\r
 314                         break;\r
 315 \r
 316                 case 'D':\r
 317                         assert(PLA > 0);\r
 318                         cEdgeType = TBD(PLA, PLB);\r
 319                         --PLA;\r
 320                         break;\r
 321 \r
 322                 case 'I':\r
 323                         assert(PLB > 0);\r
 324                         cEdgeType = TBI(PLA, PLB);\r
 325                         --PLB;\r
 326                         break;\r
 327                 \r
 328                 default:\r
 329                         Quit("Invalid edge %c", cEdgeType);\r
 330                         }\r
 331                 if (0 == PLA && 0 == PLB)\r
 332                         break;\r
 333                 }\r
 334         Path.Validate();\r
 335 \r
 336 //      SCORE Score = TraceBack(PA, uLengthA, PB, uLengthB, DPM_, DPD_, DPI_, Path);\r
 337 \r
 338 #if     TRACE\r
 339         SCORE scorePath = FastScorePath2(PA, uLengthA, PB, uLengthB, Path);\r
 340         Path.LogMe();\r
 341         Log("Score = %s Path = %s\n", LocalScoreToStr(BestScore), LocalScoreToStr(scorePath));\r
 342 #endif\r
 343 \r
 344         if (g_bKeepSimpleDP)\r
 345                 {\r
 346                 g_DPM = DPM_;\r
 347                 g_DPD = DPD_;\r
 348                 g_DPI = DPI_;\r
 349 \r
 350                 g_TBM = TBM_;\r
 351                 g_TBD = TBD_;\r
 352                 g_TBI = TBI_;\r
 353                 }\r
 354         else\r
 355                 {\r
 356                 delete[] DPM_;\r
 357                 delete[] DPD_;\r
 358                 delete[] DPI_;\r
 359 \r
 360                 delete[] TBM_;\r
 361                 delete[] TBD_;\r
 362                 delete[] TBI_;\r
 363                 }\r
 364 \r
 365         return BestScore;\r
 366         }\r
 367 \r
 368 #endif // SINLGLE_AFFINE\r