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

   1 #include "muscle.h"\r
   2 #include "profile.h"\r
   3 #include "pwpath.h"\r
   4 #include "seq.h"\r
   5 \r
   6 extern SCOREMATRIX VTML_SP;\r
   7 \r
   8 // #define SUBST(i, j)  Subst(seqA, seqB, i, j)\r
   9 #define SUBST(i, j)             MxRowA[i][seqB.GetLetter(j)]\r
  10 \r
  11 static SCORE Subst(const Seq &seqA, const Seq &seqB, unsigned i, unsigned j)\r
  12         {\r
  13         assert(i < seqA.Length());\r
  14         assert(j < seqB.Length());\r
  15 \r
  16         unsigned uLetterA = seqA.GetLetter(i);\r
  17         unsigned uLetterB = seqB.GetLetter(j);\r
  18         return VTML_SP[uLetterA][uLetterB] + g_scoreCenter;\r
  19         }\r
  20 \r
  21 struct DP_MEMORY\r
  22         {\r
  23         unsigned uLength;\r
  24         SCORE *MPrev;\r
  25         SCORE *MCurr;\r
  26         SCORE *MWork;\r
  27         SCORE *DPrev;\r
  28         SCORE *DCurr;\r
  29         SCORE *DWork;\r
  30         SCORE **MxRowA;\r
  31         unsigned *LettersB;\r
  32         unsigned *uDeletePos;\r
  33         int **TraceBack;\r
  34         };\r
  35 \r
  36 static struct DP_MEMORY DPM;\r
  37 \r
  38 static void AllocDPMem(unsigned uLengthA, unsigned uLengthB)\r
  39         {\r
  40 // Max prefix length\r
  41         unsigned uLength = (uLengthA > uLengthB ? uLengthA : uLengthB) + 1;\r
  42         if (uLength < DPM.uLength)\r
  43                 return;\r
  44 \r
  45 // Add 256 to allow for future expansion and\r
  46 // round up to next multiple of 32.\r
  47         uLength += 256;\r
  48         uLength += 32 - uLength%32;\r
  49 \r
  50         const unsigned uOldLength = DPM.uLength;\r
  51         if (uOldLength > 0)\r
  52                 {\r
  53                 for (unsigned i = 0; i < uOldLength; ++i)\r
  54                         delete[] DPM.TraceBack[i];\r
  55 \r
  56                 delete[] DPM.MPrev;\r
  57                 delete[] DPM.MCurr;\r
  58                 delete[] DPM.MWork;\r
  59                 delete[] DPM.DPrev;\r
  60                 delete[] DPM.DCurr;\r
  61                 delete[] DPM.DWork;\r
  62                 delete[] DPM.MxRowA;\r
  63                 delete[] DPM.LettersB;\r
  64                 delete[] DPM.uDeletePos;\r
  65                 delete[] DPM.TraceBack;\r
  66                 }\r
  67 \r
  68         DPM.uLength = uLength;\r
  69 \r
  70         DPM.MPrev = new SCORE[uLength];\r
  71         DPM.MCurr = new SCORE[uLength];\r
  72         DPM.MWork = new SCORE[uLength];\r
  73 \r
  74         DPM.DPrev = new SCORE[uLength];\r
  75         DPM.DCurr = new SCORE[uLength];\r
  76         DPM.DWork = new SCORE[uLength];\r
  77         DPM.MxRowA = new SCORE *[uLength];\r
  78         DPM.LettersB = new unsigned[uLength];\r
  79         DPM.uDeletePos = new unsigned[uLength];\r
  80 \r
  81         DPM.TraceBack = new int*[uLength];\r
  82 \r
  83         for (unsigned i = 0; i < uLength; ++i)\r
  84                 DPM.TraceBack[i] = new int[uLength];\r
  85         }\r
  86 \r
  87 static void RowFromSeq(const Seq &s, SCORE *Row[])\r
  88         {\r
  89         const unsigned uLength = s.Length();\r
  90         for (unsigned i = 0; i < uLength; ++i)\r
  91                 {\r
  92                 char c = s.GetChar(i);\r
  93                 unsigned uLetter = CharToLetter(c);\r
  94                 if (uLetter < 20)\r
  95                         Row[i] = VTML_SP[uLetter];\r
  96                 else\r
  97                         Row[i] = VTML_SP[AX_X];\r
  98                 }\r
  99         }\r
 100 \r
 101 static void LettersFromSeq(const Seq &s, unsigned Letters[])\r
 102         {\r
 103         const unsigned uLength = s.Length();\r
 104         for (unsigned i = 0; i < uLength; ++i)\r
 105                 {\r
 106                 char c = s.GetChar(i);\r
 107                 unsigned uLetter = CharToLetter(c);\r
 108                 if (uLetter < 20)\r
 109                         Letters[i] = uLetter;\r
 110                 else\r
 111                         Letters[i] = AX_X;\r
 112                 }\r
 113         }\r
 114 \r
 115 SCORE GlobalAlignSS(const Seq &seqA, const Seq &seqB, PWPath &Path)\r
 116         {\r
 117         const unsigned uLengthA = seqA.Length();\r
 118         const unsigned uLengthB = seqB.Length();\r
 119         const unsigned uPrefixCountA = uLengthA + 1;\r
 120         const unsigned uPrefixCountB = uLengthB + 1;\r
 121 \r
 122         AllocDPMem(uLengthA, uLengthB);\r
 123 \r
 124         SCORE *MPrev = DPM.MPrev;\r
 125         SCORE *MCurr = DPM.MCurr;\r
 126         SCORE *MWork = DPM.MWork;\r
 127 \r
 128         SCORE *DPrev = DPM.DPrev;\r
 129         SCORE *DCurr = DPM.DCurr;\r
 130         SCORE *DWork = DPM.DWork;\r
 131         SCORE **MxRowA = DPM.MxRowA;\r
 132         unsigned *LettersB = DPM.LettersB;\r
 133 \r
 134         RowFromSeq(seqA, MxRowA);\r
 135         LettersFromSeq(seqB, LettersB);\r
 136 \r
 137         unsigned *uDeletePos = DPM.uDeletePos;\r
 138 \r
 139         int **TraceBack = DPM.TraceBack;\r
 140 \r
 141 #if     DEBUG\r
 142         for (unsigned i = 0; i < uPrefixCountA; ++i)\r
 143                 memset(TraceBack[i], 0, uPrefixCountB*sizeof(int));\r
 144 #endif\r
 145 \r
 146 // Special case for i=0\r
 147         TraceBack[0][0] = 0;\r
 148         MPrev[0] = MxRowA[0][LettersB[0]];\r
 149 \r
 150 // D(0,0) is -infinity (requires I->D).\r
 151         DPrev[0] = MINUS_INFINITY;\r
 152 \r
 153         for (unsigned j = 1; j < uLengthB; ++j)\r
 154                 {\r
 155                 unsigned uLetterB = LettersB[j];\r
 156 \r
 157         // Only way to get M(0, j) looks like this:\r
 158         //              A       ----X\r
 159         //              B       XXXXX\r
 160         //                      0   j\r
 161         // So gap-open at j=0, gap-close at j-1.\r
 162                 MPrev[j] = MxRowA[0][uLetterB] + g_scoreGapOpen/2; // term gaps half\r
 163                 TraceBack[0][j] = -(int) j;\r
 164 \r
 165         // Assume no D->I transitions, then can't be a delete if only\r
 166         // one letter from A.\r
 167                 DPrev[j] = MINUS_INFINITY;\r
 168                 }\r
 169 \r
 170         SCORE IPrev_j_1;\r
 171         for (unsigned i = 1; i < uLengthA; ++i)\r
 172                 {\r
 173                 SCORE *ptrMCurr_j = MCurr;\r
 174                 memset(ptrMCurr_j, 0, uLengthB*sizeof(SCORE));\r
 175 \r
 176                 const SCORE *RowA = MxRowA[i];\r
 177                 const SCORE *ptrRowA = MxRowA[i];\r
 178                 const SCORE *ptrMCurrEnd = ptrMCurr_j + uLengthB;\r
 179                 unsigned *ptrLettersB = LettersB;\r
 180                 for (; ptrMCurr_j != ptrMCurrEnd; ++ptrMCurr_j)\r
 181                         {\r
 182                         *ptrMCurr_j = RowA[*ptrLettersB];\r
 183                         ++ptrLettersB;\r
 184                         }\r
 185 \r
 186                 unsigned *ptrDeletePos = uDeletePos;\r
 187 \r
 188         // Special case for j=0\r
 189         // Only way to get M(i, 0) looks like this:\r
 190         //                      0   i\r
 191         //              A       XXXXX\r
 192         //              B       ----X\r
 193         // So gap-open at i=0, gap-close at i-1.\r
 194                 ptrMCurr_j = MCurr;\r
 195                 assert(ptrMCurr_j == &(MCurr[0]));\r
 196                 *ptrMCurr_j += g_scoreGapOpen/2;        // term gaps half\r
 197 \r
 198                 ++ptrMCurr_j;\r
 199 \r
 200                 int *ptrTraceBack_ij = TraceBack[i];\r
 201                 *ptrTraceBack_ij++ = (int) i;\r
 202 \r
 203                 SCORE *ptrMPrev_j = MPrev;\r
 204                 SCORE *ptrDPrev = DPrev;\r
 205                 SCORE d = *ptrDPrev;\r
 206                 SCORE DNew = *ptrMPrev_j + g_scoreGapOpen;\r
 207                 if (DNew > d)\r
 208                         {\r
 209                         d = DNew;\r
 210                         *ptrDeletePos = i;\r
 211                         }\r
 212 \r
 213                 SCORE *ptrDCurr = DCurr;\r
 214 \r
 215                 assert(ptrDCurr == &(DCurr[0]));\r
 216                 *ptrDCurr = d;\r
 217 \r
 218         // Can't have an insert if no letters from B\r
 219                 IPrev_j_1 = MINUS_INFINITY;\r
 220 \r
 221                 unsigned uInsertPos;\r
 222 \r
 223                 for (unsigned j = 1; j < uLengthB; ++j)\r
 224                         {\r
 225                 // Here, MPrev_j is preserved from previous\r
 226                 // iteration so with current i,j is M[i-1][j-1]\r
 227                         SCORE MPrev_j = *ptrMPrev_j;\r
 228                         SCORE INew = MPrev_j + g_scoreGapOpen;\r
 229                         if (INew > IPrev_j_1)\r
 230                                 {\r
 231                                 IPrev_j_1 = INew;\r
 232                                 uInsertPos = j;\r
 233                                 }\r
 234 \r
 235                         SCORE scoreMax = MPrev_j;\r
 236 \r
 237                         assert(ptrDPrev == &(DPrev[j-1]));\r
 238                         SCORE scoreD = *ptrDPrev++;\r
 239                         if (scoreD > scoreMax)\r
 240                                 {\r
 241                                 scoreMax = scoreD;\r
 242                                 assert(ptrDeletePos == &(uDeletePos[j-1]));\r
 243                                 *ptrTraceBack_ij = (int) i - (int) *ptrDeletePos;\r
 244                                 assert(*ptrTraceBack_ij > 0);\r
 245                                 }\r
 246                         ++ptrDeletePos;\r
 247 \r
 248                         SCORE scoreI = IPrev_j_1;\r
 249                         if (scoreI > scoreMax)\r
 250                                 {\r
 251                                 scoreMax = scoreI;\r
 252                                 *ptrTraceBack_ij = (int) uInsertPos - (int) j;\r
 253                                 assert(*ptrTraceBack_ij < 0);\r
 254                                 }\r
 255 \r
 256                         *ptrMCurr_j += scoreMax;\r
 257                         assert(ptrMCurr_j == &(MCurr[j]));\r
 258                         ++ptrMCurr_j;\r
 259 \r
 260                         MPrev_j = *(++ptrMPrev_j);\r
 261                         assert(ptrDPrev == &(DPrev[j]));\r
 262                         SCORE d = *ptrDPrev;\r
 263                         SCORE DNew = MPrev_j + g_scoreGapOpen;\r
 264                         if (DNew > d)\r
 265                                 {\r
 266                                 d = DNew;\r
 267                                 assert(ptrDeletePos == &uDeletePos[j]);\r
 268                                 *ptrDeletePos = i;\r
 269                                 }\r
 270                         assert(ptrDCurr + 1 == &(DCurr[j]));\r
 271                         *(++ptrDCurr) = d;\r
 272 \r
 273                         ++ptrTraceBack_ij;\r
 274                         }\r
 275 \r
 276                 Rotate(MPrev, MCurr, MWork);\r
 277                 Rotate(DPrev, DCurr, DWork);\r
 278                 }\r
 279 \r
 280 // Special case for i=uLengthA\r
 281         SCORE IPrev = MINUS_INFINITY;\r
 282 \r
 283         unsigned uInsertPos;\r
 284 \r
 285         for (unsigned j = 1; j < uLengthB; ++j)\r
 286                 {\r
 287                 SCORE INew = MPrev[j-1];\r
 288                 if (INew > IPrev)\r
 289                         {\r
 290                         uInsertPos = j;\r
 291                         IPrev = INew;\r
 292                         }\r
 293                 }\r
 294 \r
 295 // Special case for i=uLengthA, j=uLengthB\r
 296         SCORE scoreMax = MPrev[uLengthB-1];\r
 297         int iTraceBack = 0;\r
 298 \r
 299         SCORE scoreD = DPrev[uLengthB-1] - g_scoreGapOpen/2;    // term gaps half\r
 300         if (scoreD > scoreMax)\r
 301                 {\r
 302                 scoreMax = scoreD;\r
 303                 iTraceBack = (int) uLengthA - (int) uDeletePos[uLengthB-1];\r
 304                 }\r
 305 \r
 306         SCORE scoreI = IPrev - g_scoreGapOpen/2;\r
 307         if (scoreI > scoreMax)\r
 308                 {\r
 309                 scoreMax = scoreI;\r
 310                 iTraceBack = (int) uInsertPos - (int) uLengthB;\r
 311                 }\r
 312 \r
 313         TraceBack[uLengthA][uLengthB] = iTraceBack;\r
 314 \r
 315         TraceBackToPath(TraceBack, uLengthA, uLengthB, Path);\r
 316 \r
 317         return scoreMax;\r
 318         }\r