binaries/src/muscle/traceback.cpp

   1 #include "muscle.h"\r
   2 #include "profile.h"\r
   3 #include "pwpath.h"\r
   4 #include <math.h>\r
   5 \r
   6 #define TRACE   0\r
   7 \r
   8 #define EQ(a, b)        (fabs(a-b) < 0.1)\r
   9 \r
  10 SCORE TraceBack(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,\r
  11   unsigned uLengthB, const SCORE *DPM_, const SCORE *DPD_, const SCORE *DPI_,\r
  12   PWPath &Path)\r
  13         {\r
  14 #if     TRACE\r
  15         Log("\n");\r
  16         Log("TraceBack LengthA=%u LengthB=%u\n", uLengthA, uLengthB);\r
  17 #endif\r
  18         assert(uLengthB > 0 && uLengthA > 0);\r
  19 \r
  20         const unsigned uPrefixCountA = uLengthA + 1;\r
  21         const unsigned uPrefixCountB = uLengthB + 1;\r
  22 \r
  23         Path.Clear();\r
  24 \r
  25         unsigned uPrefixLengthA = uLengthA;\r
  26         unsigned uPrefixLengthB = uLengthB;\r
  27 \r
  28         const SCORE scoreM = DPM(uPrefixLengthA, uPrefixLengthB);\r
  29         SCORE scoreD = DPD(uPrefixLengthA, uPrefixLengthB);\r
  30         SCORE scoreI = DPI(uPrefixLengthA, uPrefixLengthB);\r
  31 \r
  32         const ProfPos &LastPPA = PA[uLengthA - 1];\r
  33         const ProfPos &LastPPB = PB[uLengthB - 1];\r
  34 \r
  35         scoreD += LastPPA.m_scoreGapClose;\r
  36         scoreI += LastPPB.m_scoreGapClose;\r
  37 \r
  38         char cEdgeType = cInsane;\r
  39         SCORE scoreMax;\r
  40         if (scoreM >= scoreD && scoreM >= scoreI)\r
  41                 {\r
  42                 scoreMax = scoreM;\r
  43                 cEdgeType = 'M';\r
  44                 }\r
  45         else if (scoreD >= scoreM && scoreD >= scoreI)\r
  46                 {\r
  47                 scoreMax = scoreD;\r
  48                 cEdgeType = 'D';\r
  49                 }\r
  50         else\r
  51                 {\r
  52                 assert(scoreI >= scoreM && scoreI >= scoreD);\r
  53                 scoreMax = scoreI;\r
  54                 cEdgeType = 'I';\r
  55                 }\r
  56 \r
  57         for (;;)\r
  58                 {\r
  59                 if ('S' == cEdgeType)\r
  60                         break;\r
  61 \r
  62                 PWEdge Edge;\r
  63                 Edge.cType = cEdgeType;\r
  64                 Edge.uPrefixLengthA = uPrefixLengthA;\r
  65                 Edge.uPrefixLengthB = uPrefixLengthB;\r
  66                 Path.PrependEdge(Edge);\r
  67 \r
  68                 char cPrevEdgeType;\r
  69                 unsigned uPrevPrefixLengthA = uPrefixLengthA;\r
  70                 unsigned uPrevPrefixLengthB = uPrefixLengthB;\r
  71 \r
  72                 switch (cEdgeType)\r
  73                         {\r
  74                 case 'M':\r
  75                         {\r
  76                         assert(uPrefixLengthA > 0);\r
  77                         assert(uPrefixLengthB > 0);\r
  78                         const ProfPos &PPA = PA[uPrefixLengthA - 1];\r
  79                         const ProfPos &PPB = PB[uPrefixLengthB - 1];\r
  80 \r
  81                         const SCORE Score = DPM(uPrefixLengthA, uPrefixLengthB);\r
  82                         const SCORE scoreMatch = ScoreProfPos2(PPA, PPB);\r
  83 \r
  84                         SCORE scoreSM;\r
  85                         if (1 == uPrefixLengthA && 1 == uPrefixLengthB)\r
  86                                 scoreSM = scoreMatch;\r
  87                         else\r
  88                                 scoreSM = MINUS_INFINITY;\r
  89 \r
  90                         SCORE scoreMM = MINUS_INFINITY;\r
  91                         SCORE scoreDM = MINUS_INFINITY;\r
  92                         SCORE scoreIM = MINUS_INFINITY;\r
  93                         if (uPrefixLengthA > 1 && uPrefixLengthB > 1)\r
  94                                 scoreMM = DPM(uPrefixLengthA-1, uPrefixLengthB-1) + scoreMatch;\r
  95                         if (uPrefixLengthA > 1)\r
  96                                 {\r
  97                                 SCORE scoreTransDM = PA[uPrefixLengthA-2].m_scoreGapClose;\r
  98                                 scoreDM = DPD(uPrefixLengthA-1, uPrefixLengthB-1) + scoreTransDM + scoreMatch;\r
  99                                 }\r
 100                         if (uPrefixLengthB > 1)\r
 101                                 {\r
 102                                 SCORE scoreTransIM = PB[uPrefixLengthB-2].m_scoreGapClose;\r
 103                                 scoreIM = DPI(uPrefixLengthA-1, uPrefixLengthB-1) + scoreTransIM + scoreMatch;\r
 104                                 }\r
 105 \r
 106                         if (EQ(scoreMM, Score))\r
 107                                 cPrevEdgeType = 'M';\r
 108                         else if (EQ(scoreDM, Score))\r
 109                                 cPrevEdgeType = 'D';\r
 110                         else if (EQ(scoreIM, Score))\r
 111                                 cPrevEdgeType = 'I';\r
 112                         else if (EQ(scoreSM, Score))\r
 113                                 cPrevEdgeType = 'S';\r
 114                         else\r
 115                                 Quit("TraceBack: failed to match M score=%g M=%g D=%g I=%g S=%g",\r
 116                                   Score, scoreMM, scoreDM, scoreIM, scoreSM);\r
 117 \r
 118                         --uPrevPrefixLengthA;\r
 119                         --uPrevPrefixLengthB;\r
 120                         break;\r
 121                         }\r
 122 \r
 123                 case 'D':\r
 124                         {\r
 125                         assert(uPrefixLengthA > 0);\r
 126                         const SCORE Score = DPD(uPrefixLengthA, uPrefixLengthB);\r
 127 \r
 128                         SCORE scoreMD = MINUS_INFINITY;\r
 129                         SCORE scoreDD = MINUS_INFINITY;\r
 130                         SCORE scoreSD = MINUS_INFINITY;\r
 131                         if (uPrefixLengthB == 0)\r
 132                                 {\r
 133                                 if (uPrefixLengthA == 1)\r
 134                                         scoreSD = PA[0].m_scoreGapOpen;\r
 135                                 else\r
 136                                         scoreSD = DPD(uPrefixLengthA - 1, 0);\r
 137                                 }\r
 138                         if (uPrefixLengthA > 1)\r
 139                                 {\r
 140                                 const ProfPos &PPA = PA[uPrefixLengthA - 1];\r
 141                                 SCORE scoreTransMD = PPA.m_scoreGapOpen;\r
 142                                 scoreMD = DPM(uPrefixLengthA-1, uPrefixLengthB) + scoreTransMD;\r
 143                                 scoreDD = DPD(uPrefixLengthA-1, uPrefixLengthB);\r
 144                                 }\r
 145 \r
 146                         if (EQ(Score, scoreMD))\r
 147                                 cPrevEdgeType = 'M';\r
 148                         else if (EQ(Score, scoreDD))\r
 149                                 cPrevEdgeType = 'D';\r
 150                         else if (EQ(Score, scoreSD))\r
 151                                 cPrevEdgeType = 'S';\r
 152                         else\r
 153                                 Quit("TraceBack: failed to match D");\r
 154 \r
 155                         --uPrevPrefixLengthA;\r
 156                         break;\r
 157                         }\r
 158 \r
 159                 case 'I':\r
 160                         {\r
 161                         assert(uPrefixLengthB > 0);\r
 162                         const SCORE Score = DPI(uPrefixLengthA, uPrefixLengthB);\r
 163 \r
 164                         SCORE scoreMI = MINUS_INFINITY;\r
 165                         SCORE scoreII = MINUS_INFINITY;\r
 166                         SCORE scoreSI = MINUS_INFINITY;\r
 167                         if (uPrefixLengthA == 0)\r
 168                                 {\r
 169                                 if (uPrefixLengthB == 1)\r
 170                                         scoreSI = PB[0].m_scoreGapOpen;\r
 171                                 else\r
 172                                         scoreSI = DPI(0, uPrefixLengthB - 1);\r
 173                                 }\r
 174                         if (uPrefixLengthB > 1)\r
 175                                 {\r
 176                                 const ProfPos &PPB = PB[uPrefixLengthB - 1];\r
 177                                 SCORE scoreTransMI = PPB.m_scoreGapOpen;\r
 178                                 scoreMI = DPM(uPrefixLengthA, uPrefixLengthB-1) + scoreTransMI;\r
 179                                 scoreII = DPI(uPrefixLengthA, uPrefixLengthB-1);\r
 180                                 }\r
 181 \r
 182                         if (EQ(Score, scoreMI))\r
 183                                 cPrevEdgeType = 'M';\r
 184                         else if (EQ(Score, scoreII))\r
 185                                 cPrevEdgeType = 'I';\r
 186                         else if (EQ(Score, scoreSI))\r
 187                                 cPrevEdgeType = 'S';\r
 188                         else\r
 189                                 Quit("TraceBack: failed to match I");\r
 190 \r
 191                         --uPrevPrefixLengthB;\r
 192                         break;\r
 193                         }\r
 194 \r
 195                 default:\r
 196                         assert(false);\r
 197                         }\r
 198 #if     TRACE\r
 199                 Log("Edge %c%c%u.%u", cPrevEdgeType, cEdgeType, uPrefixLengthA, uPrefixLengthB);\r
 200                 Log("\n");\r
 201 #endif\r
 202                 cEdgeType = cPrevEdgeType;\r
 203                 uPrefixLengthA = uPrevPrefixLengthA;\r
 204                 uPrefixLengthB = uPrevPrefixLengthB;\r
 205                 }\r
 206 \r
 207         return scoreMax;\r
 208         }\r