binaries/src/muscle/sw.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 // Textbook Smith-Waterman affine gap implementation.\r
   8 \r
   9 #define TRACE   0\r
  10 \r
  11 static const char *LocalScoreToStr(SCORE s)\r
  12         {\r
  13         static char str[16];\r
  14         if (MINUS_INFINITY == s)\r
  15                 return "     *";\r
  16         sprintf(str, "%6.2f", s);\r
  17         return str;\r
  18         }\r
  19 \r
  20 static void ListDP(const SCORE *DPM_, const ProfPos *PA, const ProfPos *PB,\r
  21   unsigned uPrefixCountA, unsigned uPrefixCountB)\r
  22         {\r
  23         Log("        ");\r
  24         for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
  25                 {\r
  26                 char c = ' ';\r
  27                 if (uPrefixLengthB > 0)\r
  28                         c = ConsensusChar(PB[uPrefixLengthB - 1]);\r
  29                 Log(" %4u:%c", uPrefixLengthB, c);\r
  30                 }\r
  31         Log("\n");\r
  32         for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)\r
  33                 {\r
  34                 char c = ' ';\r
  35                 if (uPrefixLengthA > 0)\r
  36                         c = ConsensusChar(PA[uPrefixLengthA - 1]);\r
  37                 Log("%4u:%c  ", uPrefixLengthA, c);\r
  38                 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
  39                         Log(" %s", LocalScoreToStr(DPM(uPrefixLengthA, uPrefixLengthB)));\r
  40                 Log("\n");\r
  41                 }\r
  42         }\r
  43 \r
  44 SCORE SW(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,\r
  45   unsigned uLengthB, PWPath &Path)\r
  46         {\r
  47         assert(uLengthB > 0 && uLengthA > 0);\r
  48 \r
  49         const unsigned uPrefixCountA = uLengthA + 1;\r
  50         const unsigned uPrefixCountB = uLengthB + 1;\r
  51 \r
  52 // Allocate DP matrices\r
  53         const size_t LM = uPrefixCountA*uPrefixCountB;\r
  54         SCORE *DPM_ = new SCORE[LM];\r
  55         SCORE *DPD_ = new SCORE[LM];\r
  56         SCORE *DPI_ = new SCORE[LM];\r
  57 \r
  58         DPM(0, 0) = 0;\r
  59         DPD(0, 0) = MINUS_INFINITY;\r
  60         DPI(0, 0) = MINUS_INFINITY;\r
  61 \r
  62         DPM(1, 0) = MINUS_INFINITY;\r
  63         DPD(1, 0) = MINUS_INFINITY;\r
  64         DPI(1, 0) = MINUS_INFINITY;\r
  65 \r
  66         DPM(0, 1) = MINUS_INFINITY;\r
  67         DPD(0, 1) = MINUS_INFINITY;\r
  68         DPI(0, 1) = MINUS_INFINITY;\r
  69 \r
  70 // Empty prefix of B is special case\r
  71         for (unsigned uPrefixLengthA = 2; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)\r
  72                 {\r
  73         // M=LetterA+LetterB, impossible with empty prefix\r
  74                 DPM(uPrefixLengthA, 0) = MINUS_INFINITY;\r
  75 \r
  76         // D=LetterA+GapB, never optimal in local alignment with gap penalties\r
  77                 DPD(uPrefixLengthA, 0) = MINUS_INFINITY;\r
  78 \r
  79         // I=GapA+LetterB, impossible with empty prefix\r
  80                 DPI(uPrefixLengthA, 0) = MINUS_INFINITY;\r
  81                 }\r
  82 \r
  83 // Empty prefix of A is special case\r
  84         for (unsigned uPrefixLengthB = 2; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
  85                 {\r
  86         // M=LetterA+LetterB, impossible with empty prefix\r
  87                 DPM(0, uPrefixLengthB) = MINUS_INFINITY;\r
  88 \r
  89         // D=LetterA+GapB, impossible with empty prefix\r
  90                 DPD(0, uPrefixLengthB) = MINUS_INFINITY;\r
  91 \r
  92         // I=GapA+LetterB, never optimal in local alignment with gap penalties\r
  93                 DPI(0, uPrefixLengthB) = MINUS_INFINITY;\r
  94                 }\r
  95 \r
  96         SCORE scoreMax = MINUS_INFINITY;\r
  97         unsigned uPrefixLengthAMax = uInsane;\r
  98         unsigned uPrefixLengthBMax = uInsane;\r
  99 \r
 100 // ============\r
 101 // Main DP loop\r
 102 // ============\r
 103         SCORE scoreGapCloseB = MINUS_INFINITY;\r
 104         for (unsigned uPrefixLengthB = 1; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
 105                 {\r
 106                 const ProfPos &PPB = PB[uPrefixLengthB - 1];\r
 107 \r
 108                 SCORE scoreGapCloseA = MINUS_INFINITY;\r
 109                 for (unsigned uPrefixLengthA = 1; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)\r
 110                         {\r
 111                         const ProfPos &PPA = PA[uPrefixLengthA - 1];\r
 112 \r
 113                         {\r
 114                 // Match M=LetterA+LetterB\r
 115                         SCORE scoreLL = ScoreProfPos2(PPA, PPB);\r
 116 \r
 117                         SCORE scoreMM = DPM(uPrefixLengthA-1, uPrefixLengthB-1);\r
 118                         SCORE scoreDM = DPD(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseA;\r
 119                         SCORE scoreIM = DPI(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseB;\r
 120 \r
 121                         SCORE scoreBest;\r
 122                         if (scoreMM >= scoreDM && scoreMM >= scoreIM)\r
 123                                 scoreBest = scoreMM;\r
 124                         else if (scoreDM >= scoreMM && scoreDM >= scoreIM)\r
 125                                 scoreBest = scoreDM;\r
 126                         else \r
 127                                 {\r
 128                                 assert(scoreIM >= scoreMM && scoreIM >= scoreDM);\r
 129                                 scoreBest = scoreIM;\r
 130                                 }\r
 131                         if (scoreBest < 0)\r
 132                                 scoreBest = 0;\r
 133                         scoreBest += scoreLL;\r
 134                         if (scoreBest > scoreMax)\r
 135                                 {\r
 136                                 scoreMax = scoreBest;\r
 137                                 uPrefixLengthAMax = uPrefixLengthA;\r
 138                                 uPrefixLengthBMax = uPrefixLengthB;\r
 139                                 }\r
 140                         DPM(uPrefixLengthA, uPrefixLengthB) = scoreBest;\r
 141                         }\r
 142 \r
 143                         {\r
 144                 // Delete D=LetterA+GapB\r
 145                         SCORE scoreMD = DPM(uPrefixLengthA-1, uPrefixLengthB) +\r
 146                           PA[uPrefixLengthA-1].m_scoreGapOpen;\r
 147                         SCORE scoreDD = DPD(uPrefixLengthA-1, uPrefixLengthB);\r
 148 \r
 149                         SCORE scoreBest;\r
 150                         if (scoreMD >= scoreDD)\r
 151                                 scoreBest = scoreMD;\r
 152                         else\r
 153                                 {\r
 154                                 assert(scoreDD >= scoreMD);\r
 155                                 scoreBest = scoreDD;\r
 156                                 }\r
 157                         DPD(uPrefixLengthA, uPrefixLengthB) = scoreBest;\r
 158                         }\r
 159 \r
 160                 // Insert I=GapA+LetterB\r
 161                         {\r
 162                         SCORE scoreMI = DPM(uPrefixLengthA, uPrefixLengthB-1) +\r
 163                           PB[uPrefixLengthB - 1].m_scoreGapOpen;\r
 164                         SCORE scoreII = DPI(uPrefixLengthA, uPrefixLengthB-1);\r
 165 \r
 166                         SCORE scoreBest;\r
 167                         if (scoreMI >= scoreII)\r
 168                                 scoreBest = scoreMI;\r
 169                         else \r
 170                                 {\r
 171                                 assert(scoreII > scoreMI);\r
 172                                 scoreBest = scoreII;\r
 173                                 }\r
 174                         DPI(uPrefixLengthA, uPrefixLengthB) = scoreBest;\r
 175                         }\r
 176 \r
 177                         scoreGapCloseA = PPA.m_scoreGapClose;\r
 178                         }\r
 179                 scoreGapCloseB = PPB.m_scoreGapClose;\r
 180                 }\r
 181 \r
 182 #if TRACE\r
 183         Log("DPM:\n");\r
 184         ListDP(DPM_, PA, PB, uPrefixLengthA, uPrefixLengthB);\r
 185         Log("DPD:\n");\r
 186         ListDP(DPD_, PA, PB, uPrefixLengthA, uPrefixLengthB);\r
 187         Log("DPI:\n");\r
 188         ListDP(DPI_, PA, PB, uPrefixLengthA, uPrefixLengthB);\r
 189 #endif\r
 190 \r
 191         assert(scoreMax == DPM(uPrefixLengthAMax, uPrefixLengthBMax));\r
 192         TraceBackSW(PA, uLengthA, PB, uLengthB, DPM_, DPD_, DPI_, \r
 193           uPrefixLengthAMax, uPrefixLengthBMax, Path);\r
 194 \r
 195 #if     TRACE\r
 196         SCORE scorePath = FastScorePath2(PA, uLengthA, PB, uLengthB, Path);\r
 197         Path.LogMe();\r
 198         Log("Score = %s Path = %s\n", LocalScoreToStr(scoreMax), LocalScoreToStr(scorePath));\r
 199 #endif\r
 200 \r
 201         delete[] DPM_;\r
 202         delete[] DPD_;\r
 203         delete[] DPI_;\r
 204 \r
 205         return scoreMax;\r
 206         }\r