| 1 | #include "muscle.h" |
|---|
| 2 | #include "profile.h" |
|---|
| 3 | #include "pwpath.h" |
|---|
| 4 | |
|---|
| 5 | struct DP_MEMORY |
|---|
| 6 | { |
|---|
| 7 | unsigned uLength; |
|---|
| 8 | SCORE *GapOpenA; |
|---|
| 9 | SCORE *GapOpenB; |
|---|
| 10 | SCORE *GapCloseA; |
|---|
| 11 | SCORE *GapCloseB; |
|---|
| 12 | SCORE *MPrev; |
|---|
| 13 | SCORE *MCurr; |
|---|
| 14 | SCORE *MWork; |
|---|
| 15 | SCORE *DPrev; |
|---|
| 16 | SCORE *DCurr; |
|---|
| 17 | SCORE *DWork; |
|---|
| 18 | SCORE **ScoreMxB; |
|---|
| 19 | unsigned **SortOrderA; |
|---|
| 20 | unsigned *uDeletePos; |
|---|
| 21 | FCOUNT **FreqsA; |
|---|
| 22 | int **TraceBack; |
|---|
| 23 | }; |
|---|
| 24 | |
|---|
| 25 | static struct DP_MEMORY DPM; |
|---|
| 26 | |
|---|
| 27 | void FreeDPMemSPN() |
|---|
| 28 | { |
|---|
| 29 | const unsigned uOldLength = DPM.uLength; |
|---|
| 30 | if (0 == uOldLength) |
|---|
| 31 | return; |
|---|
| 32 | |
|---|
| 33 | for (unsigned i = 0; i < uOldLength; ++i) |
|---|
| 34 | { |
|---|
| 35 | delete[] DPM.TraceBack[i]; |
|---|
| 36 | delete[] DPM.FreqsA[i]; |
|---|
| 37 | delete[] DPM.SortOrderA[i]; |
|---|
| 38 | } |
|---|
| 39 | for (unsigned n = 0; n < 4; ++n) |
|---|
| 40 | delete[] DPM.ScoreMxB[n]; |
|---|
| 41 | |
|---|
| 42 | delete[] DPM.MPrev; |
|---|
| 43 | delete[] DPM.MCurr; |
|---|
| 44 | delete[] DPM.MWork; |
|---|
| 45 | delete[] DPM.DPrev; |
|---|
| 46 | delete[] DPM.DCurr; |
|---|
| 47 | delete[] DPM.DWork; |
|---|
| 48 | delete[] DPM.uDeletePos; |
|---|
| 49 | delete[] DPM.GapOpenA; |
|---|
| 50 | delete[] DPM.GapOpenB; |
|---|
| 51 | delete[] DPM.GapCloseA; |
|---|
| 52 | delete[] DPM.GapCloseB; |
|---|
| 53 | delete[] DPM.SortOrderA; |
|---|
| 54 | delete[] DPM.FreqsA; |
|---|
| 55 | delete[] DPM.ScoreMxB; |
|---|
| 56 | delete[] DPM.TraceBack; |
|---|
| 57 | } |
|---|
| 58 | |
|---|
| 59 | static void AllocDPMem(unsigned uLengthA, unsigned uLengthB) |
|---|
| 60 | { |
|---|
| 61 | // Max prefix length |
|---|
| 62 | unsigned uLength = (uLengthA > uLengthB ? uLengthA : uLengthB) + 1; |
|---|
| 63 | if (uLength < DPM.uLength) |
|---|
| 64 | return; |
|---|
| 65 | |
|---|
| 66 | // Add 256 to allow for future expansion and |
|---|
| 67 | // round up to next multiple of 32. |
|---|
| 68 | uLength += 256; |
|---|
| 69 | uLength += 32 - uLength%32; |
|---|
| 70 | |
|---|
| 71 | const unsigned uOldLength = DPM.uLength; |
|---|
| 72 | if (uOldLength > 0) |
|---|
| 73 | { |
|---|
| 74 | for (unsigned i = 0; i < uOldLength; ++i) |
|---|
| 75 | { |
|---|
| 76 | delete[] DPM.TraceBack[i]; |
|---|
| 77 | delete[] DPM.FreqsA[i]; |
|---|
| 78 | delete[] DPM.SortOrderA[i]; |
|---|
| 79 | } |
|---|
| 80 | for (unsigned n = 0; n < 4; ++n) |
|---|
| 81 | delete[] DPM.ScoreMxB[n]; |
|---|
| 82 | |
|---|
| 83 | delete[] DPM.MPrev; |
|---|
| 84 | delete[] DPM.MCurr; |
|---|
| 85 | delete[] DPM.MWork; |
|---|
| 86 | delete[] DPM.DPrev; |
|---|
| 87 | delete[] DPM.DCurr; |
|---|
| 88 | delete[] DPM.DWork; |
|---|
| 89 | delete[] DPM.uDeletePos; |
|---|
| 90 | delete[] DPM.GapOpenA; |
|---|
| 91 | delete[] DPM.GapOpenB; |
|---|
| 92 | delete[] DPM.GapCloseA; |
|---|
| 93 | delete[] DPM.GapCloseB; |
|---|
| 94 | delete[] DPM.SortOrderA; |
|---|
| 95 | delete[] DPM.FreqsA; |
|---|
| 96 | delete[] DPM.ScoreMxB; |
|---|
| 97 | delete[] DPM.TraceBack; |
|---|
| 98 | } |
|---|
| 99 | |
|---|
| 100 | DPM.uLength = uLength; |
|---|
| 101 | |
|---|
| 102 | DPM.GapOpenA = new SCORE[uLength]; |
|---|
| 103 | DPM.GapOpenB = new SCORE[uLength]; |
|---|
| 104 | DPM.GapCloseA = new SCORE[uLength]; |
|---|
| 105 | DPM.GapCloseB = new SCORE[uLength]; |
|---|
| 106 | |
|---|
| 107 | DPM.SortOrderA = new unsigned*[uLength]; |
|---|
| 108 | DPM.FreqsA = new FCOUNT*[uLength]; |
|---|
| 109 | DPM.ScoreMxB = new SCORE*[4]; |
|---|
| 110 | DPM.MPrev = new SCORE[uLength]; |
|---|
| 111 | DPM.MCurr = new SCORE[uLength]; |
|---|
| 112 | DPM.MWork = new SCORE[uLength]; |
|---|
| 113 | |
|---|
| 114 | DPM.DPrev = new SCORE[uLength]; |
|---|
| 115 | DPM.DCurr = new SCORE[uLength]; |
|---|
| 116 | DPM.DWork = new SCORE[uLength]; |
|---|
| 117 | DPM.uDeletePos = new unsigned[uLength]; |
|---|
| 118 | |
|---|
| 119 | DPM.TraceBack = new int*[uLength]; |
|---|
| 120 | |
|---|
| 121 | for (unsigned uLetter = 0; uLetter < 4; ++uLetter) |
|---|
| 122 | DPM.ScoreMxB[uLetter] = new SCORE[uLength]; |
|---|
| 123 | |
|---|
| 124 | for (unsigned i = 0; i < uLength; ++i) |
|---|
| 125 | { |
|---|
| 126 | DPM.SortOrderA[i] = new unsigned[4]; |
|---|
| 127 | DPM.FreqsA[i] = new FCOUNT[4]; |
|---|
| 128 | DPM.TraceBack[i] = new int[uLength]; |
|---|
| 129 | } |
|---|
| 130 | } |
|---|
| 131 | |
|---|
| 132 | SCORE GlobalAlignSPN(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB, |
|---|
| 133 | unsigned uLengthB, PWPath &Path) |
|---|
| 134 | { |
|---|
| 135 | if (ALPHA_DNA != g_Alpha || ALPHA_RNA == g_Alpha) |
|---|
| 136 | Quit("GlobalAlignSPN: must be nucleo"); |
|---|
| 137 | |
|---|
| 138 | const unsigned uPrefixCountA = uLengthA + 1; |
|---|
| 139 | const unsigned uPrefixCountB = uLengthB + 1; |
|---|
| 140 | |
|---|
| 141 | AllocDPMem(uLengthA, uLengthB); |
|---|
| 142 | |
|---|
| 143 | SCORE *GapOpenA = DPM.GapOpenA; |
|---|
| 144 | SCORE *GapOpenB = DPM.GapOpenB; |
|---|
| 145 | SCORE *GapCloseA = DPM.GapCloseA; |
|---|
| 146 | SCORE *GapCloseB = DPM.GapCloseB; |
|---|
| 147 | |
|---|
| 148 | unsigned **SortOrderA = DPM.SortOrderA; |
|---|
| 149 | FCOUNT **FreqsA = DPM.FreqsA; |
|---|
| 150 | SCORE **ScoreMxB = DPM.ScoreMxB; |
|---|
| 151 | SCORE *MPrev = DPM.MPrev; |
|---|
| 152 | SCORE *MCurr = DPM.MCurr; |
|---|
| 153 | SCORE *MWork = DPM.MWork; |
|---|
| 154 | |
|---|
| 155 | SCORE *DPrev = DPM.DPrev; |
|---|
| 156 | SCORE *DCurr = DPM.DCurr; |
|---|
| 157 | SCORE *DWork = DPM.DWork; |
|---|
| 158 | unsigned *uDeletePos = DPM.uDeletePos; |
|---|
| 159 | |
|---|
| 160 | int **TraceBack = DPM.TraceBack; |
|---|
| 161 | |
|---|
| 162 | for (unsigned i = 0; i < uLengthA; ++i) |
|---|
| 163 | { |
|---|
| 164 | GapOpenA[i] = PA[i].m_scoreGapOpen; |
|---|
| 165 | GapCloseA[i] = PA[i].m_scoreGapClose; |
|---|
| 166 | |
|---|
| 167 | for (unsigned uLetter = 0; uLetter < 4; ++uLetter) |
|---|
| 168 | { |
|---|
| 169 | SortOrderA[i][uLetter] = PA[i].m_uSortOrder[uLetter]; |
|---|
| 170 | FreqsA[i][uLetter] = PA[i].m_fcCounts[uLetter]; |
|---|
| 171 | } |
|---|
| 172 | } |
|---|
| 173 | |
|---|
| 174 | for (unsigned j = 0; j < uLengthB; ++j) |
|---|
| 175 | { |
|---|
| 176 | GapOpenB[j] = PB[j].m_scoreGapOpen; |
|---|
| 177 | GapCloseB[j] = PB[j].m_scoreGapClose; |
|---|
| 178 | } |
|---|
| 179 | |
|---|
| 180 | for (unsigned uLetter = 0; uLetter < 4; ++uLetter) |
|---|
| 181 | { |
|---|
| 182 | for (unsigned j = 0; j < uLengthB; ++j) |
|---|
| 183 | ScoreMxB[uLetter][j] = PB[j].m_AAScores[uLetter]; |
|---|
| 184 | } |
|---|
| 185 | |
|---|
| 186 | for (unsigned i = 0; i < uPrefixCountA; ++i) |
|---|
| 187 | memset(TraceBack[i], 0, uPrefixCountB*sizeof(int)); |
|---|
| 188 | |
|---|
| 189 | // Special case for i=0 |
|---|
| 190 | unsigned **ptrSortOrderA = SortOrderA; |
|---|
| 191 | FCOUNT **ptrFreqsA = FreqsA; |
|---|
| 192 | assert(ptrSortOrderA == &(SortOrderA[0])); |
|---|
| 193 | assert(ptrFreqsA == &(FreqsA[0])); |
|---|
| 194 | TraceBack[0][0] = 0; |
|---|
| 195 | |
|---|
| 196 | SCORE scoreSum = 0; |
|---|
| 197 | unsigned *ptrSortOrderAi = SortOrderA[0]; |
|---|
| 198 | const unsigned *ptrSortOrderAEnd = ptrSortOrderAi + 4; |
|---|
| 199 | FCOUNT *ptrFreqsAi = FreqsA[0]; |
|---|
| 200 | for (; ptrSortOrderAi != ptrSortOrderAEnd; ++ptrSortOrderAi) |
|---|
| 201 | { |
|---|
| 202 | const unsigned uLetter = *ptrSortOrderAi; |
|---|
| 203 | const FCOUNT fcLetter = ptrFreqsAi[uLetter]; |
|---|
| 204 | if (0 == fcLetter) |
|---|
| 205 | break; |
|---|
| 206 | scoreSum += fcLetter*ScoreMxB[uLetter][0]; |
|---|
| 207 | } |
|---|
| 208 | MPrev[0] = scoreSum - g_scoreCenter; |
|---|
| 209 | |
|---|
| 210 | // D(0,0) is -infinity (requires I->D). |
|---|
| 211 | DPrev[0] = MINUS_INFINITY; |
|---|
| 212 | |
|---|
| 213 | for (unsigned j = 1; j < uLengthB; ++j) |
|---|
| 214 | { |
|---|
| 215 | // Only way to get M(0, j) looks like this: |
|---|
| 216 | // A ----X |
|---|
| 217 | // B XXXXX |
|---|
| 218 | // 0 j |
|---|
| 219 | // So gap-open at j=0, gap-close at j-1. |
|---|
| 220 | SCORE scoreSum = 0; |
|---|
| 221 | unsigned *ptrSortOrderAi = SortOrderA[0]; |
|---|
| 222 | const unsigned *ptrSortOrderAEnd = ptrSortOrderAi + 4; |
|---|
| 223 | FCOUNT *ptrFreqsAi = FreqsA[0]; |
|---|
| 224 | for (; ptrSortOrderAi != ptrSortOrderAEnd; ++ptrSortOrderAi) |
|---|
| 225 | { |
|---|
| 226 | const unsigned uLetter = *ptrSortOrderAi; |
|---|
| 227 | const FCOUNT fcLetter = ptrFreqsAi[uLetter]; |
|---|
| 228 | if (0 == fcLetter) |
|---|
| 229 | break; |
|---|
| 230 | scoreSum += fcLetter*ScoreMxB[uLetter][j]; |
|---|
| 231 | } |
|---|
| 232 | MPrev[j] = scoreSum - g_scoreCenter + GapOpenB[0] + GapCloseB[j-1]; |
|---|
| 233 | TraceBack[0][j] = -(int) j; |
|---|
| 234 | |
|---|
| 235 | // Assume no D->I transitions, then can't be a delete if only |
|---|
| 236 | // one letter from A. |
|---|
| 237 | DPrev[j] = MINUS_INFINITY; |
|---|
| 238 | } |
|---|
| 239 | |
|---|
| 240 | SCORE IPrev_j_1; |
|---|
| 241 | for (unsigned i = 1; i < uLengthA; ++i) |
|---|
| 242 | { |
|---|
| 243 | ++ptrSortOrderA; |
|---|
| 244 | ++ptrFreqsA; |
|---|
| 245 | assert(ptrSortOrderA == &(SortOrderA[i])); |
|---|
| 246 | assert(ptrFreqsA == &(FreqsA[i])); |
|---|
| 247 | |
|---|
| 248 | SCORE *ptrMCurr_j = MCurr; |
|---|
| 249 | memset(ptrMCurr_j, 0, uLengthB*sizeof(SCORE)); |
|---|
| 250 | const FCOUNT *FreqsAi = *ptrFreqsA; |
|---|
| 251 | |
|---|
| 252 | const unsigned *SortOrderAi = *ptrSortOrderA; |
|---|
| 253 | const unsigned *ptrSortOrderAiEnd = SortOrderAi + 4; |
|---|
| 254 | const SCORE *ptrMCurrMax = MCurr + uLengthB; |
|---|
| 255 | for (const unsigned *ptrSortOrderAi = SortOrderAi; |
|---|
| 256 | ptrSortOrderAi != ptrSortOrderAiEnd; |
|---|
| 257 | ++ptrSortOrderAi) |
|---|
| 258 | { |
|---|
| 259 | const unsigned uLetter = *ptrSortOrderAi; |
|---|
| 260 | SCORE *NSBR_Letter = ScoreMxB[uLetter]; |
|---|
| 261 | const FCOUNT fcLetter = FreqsAi[uLetter]; |
|---|
| 262 | if (0 == fcLetter) |
|---|
| 263 | break; |
|---|
| 264 | SCORE *ptrNSBR = NSBR_Letter; |
|---|
| 265 | for (SCORE *ptrMCurr = MCurr; ptrMCurr != ptrMCurrMax; ++ptrMCurr) |
|---|
| 266 | *ptrMCurr += fcLetter*(*ptrNSBR++); |
|---|
| 267 | } |
|---|
| 268 | |
|---|
| 269 | for (unsigned j = 0; j < uLengthB; ++j) |
|---|
| 270 | MCurr[j] -= g_scoreCenter; |
|---|
| 271 | |
|---|
| 272 | ptrMCurr_j = MCurr; |
|---|
| 273 | unsigned *ptrDeletePos = uDeletePos; |
|---|
| 274 | |
|---|
| 275 | // Special case for j=0 |
|---|
| 276 | // Only way to get M(i, 0) looks like this: |
|---|
| 277 | // 0 i |
|---|
| 278 | // A XXXXX |
|---|
| 279 | // B ----X |
|---|
| 280 | // So gap-open at i=0, gap-close at i-1. |
|---|
| 281 | assert(ptrMCurr_j == &(MCurr[0])); |
|---|
| 282 | *ptrMCurr_j += GapOpenA[0] + GapCloseA[i-1]; |
|---|
| 283 | |
|---|
| 284 | ++ptrMCurr_j; |
|---|
| 285 | |
|---|
| 286 | int *ptrTraceBack_ij = TraceBack[i]; |
|---|
| 287 | *ptrTraceBack_ij++ = (int) i; |
|---|
| 288 | |
|---|
| 289 | SCORE *ptrMPrev_j = MPrev; |
|---|
| 290 | SCORE *ptrDPrev = DPrev; |
|---|
| 291 | SCORE d = *ptrDPrev; |
|---|
| 292 | SCORE DNew = *ptrMPrev_j + GapOpenA[i]; |
|---|
| 293 | if (DNew > d) |
|---|
| 294 | { |
|---|
| 295 | d = DNew; |
|---|
| 296 | *ptrDeletePos = i; |
|---|
| 297 | } |
|---|
| 298 | |
|---|
| 299 | SCORE *ptrDCurr = DCurr; |
|---|
| 300 | |
|---|
| 301 | assert(ptrDCurr == &(DCurr[0])); |
|---|
| 302 | *ptrDCurr = d; |
|---|
| 303 | |
|---|
| 304 | // Can't have an insert if no letters from B |
|---|
| 305 | IPrev_j_1 = MINUS_INFINITY; |
|---|
| 306 | |
|---|
| 307 | unsigned uInsertPos; |
|---|
| 308 | const SCORE scoreGapOpenAi = GapOpenA[i]; |
|---|
| 309 | const SCORE scoreGapCloseAi_1 = GapCloseA[i-1]; |
|---|
| 310 | |
|---|
| 311 | for (unsigned j = 1; j < uLengthB; ++j) |
|---|
| 312 | { |
|---|
| 313 | // Here, MPrev_j is preserved from previous |
|---|
| 314 | // iteration so with current i,j is M[i-1][j-1] |
|---|
| 315 | SCORE MPrev_j = *ptrMPrev_j; |
|---|
| 316 | SCORE INew = MPrev_j + GapOpenB[j]; |
|---|
| 317 | if (INew > IPrev_j_1) |
|---|
| 318 | { |
|---|
| 319 | IPrev_j_1 = INew; |
|---|
| 320 | uInsertPos = j; |
|---|
| 321 | } |
|---|
| 322 | |
|---|
| 323 | SCORE scoreMax = MPrev_j; |
|---|
| 324 | |
|---|
| 325 | assert(ptrDPrev == &(DPrev[j-1])); |
|---|
| 326 | SCORE scoreD = *ptrDPrev++ + scoreGapCloseAi_1; |
|---|
| 327 | if (scoreD > scoreMax) |
|---|
| 328 | { |
|---|
| 329 | scoreMax = scoreD; |
|---|
| 330 | assert(ptrDeletePos == &(uDeletePos[j-1])); |
|---|
| 331 | *ptrTraceBack_ij = (int) i - (int) *ptrDeletePos; |
|---|
| 332 | assert(*ptrTraceBack_ij > 0); |
|---|
| 333 | } |
|---|
| 334 | ++ptrDeletePos; |
|---|
| 335 | |
|---|
| 336 | SCORE scoreI = IPrev_j_1 + GapCloseB[j-1]; |
|---|
| 337 | if (scoreI > scoreMax) |
|---|
| 338 | { |
|---|
| 339 | scoreMax = scoreI; |
|---|
| 340 | *ptrTraceBack_ij = (int) uInsertPos - (int) j; |
|---|
| 341 | assert(*ptrTraceBack_ij < 0); |
|---|
| 342 | } |
|---|
| 343 | |
|---|
| 344 | assert(ptrSortOrderA == &(SortOrderA[i])); |
|---|
| 345 | assert(ptrFreqsA == &(FreqsA[i])); |
|---|
| 346 | |
|---|
| 347 | *ptrMCurr_j += scoreMax; |
|---|
| 348 | assert(ptrMCurr_j == &(MCurr[j])); |
|---|
| 349 | ++ptrMCurr_j; |
|---|
| 350 | |
|---|
| 351 | MPrev_j = *(++ptrMPrev_j); |
|---|
| 352 | assert(ptrDPrev == &(DPrev[j])); |
|---|
| 353 | SCORE d = *ptrDPrev; |
|---|
| 354 | SCORE DNew = MPrev_j + scoreGapOpenAi; |
|---|
| 355 | if (DNew > d) |
|---|
| 356 | { |
|---|
| 357 | d = DNew; |
|---|
| 358 | assert(ptrDeletePos == &uDeletePos[j]); |
|---|
| 359 | *ptrDeletePos = i; |
|---|
| 360 | } |
|---|
| 361 | assert(ptrDCurr + 1 == &(DCurr[j])); |
|---|
| 362 | *(++ptrDCurr) = d; |
|---|
| 363 | |
|---|
| 364 | ++ptrTraceBack_ij; |
|---|
| 365 | } |
|---|
| 366 | |
|---|
| 367 | Rotate(MPrev, MCurr, MWork); |
|---|
| 368 | Rotate(DPrev, DCurr, DWork); |
|---|
| 369 | } |
|---|
| 370 | |
|---|
| 371 | // Special case for i=uLengthA |
|---|
| 372 | SCORE IPrev = MINUS_INFINITY; |
|---|
| 373 | |
|---|
| 374 | unsigned uInsertPos; |
|---|
| 375 | |
|---|
| 376 | for (unsigned j = 1; j < uLengthB; ++j) |
|---|
| 377 | { |
|---|
| 378 | SCORE INew = MPrev[j-1] + GapOpenB[j]; |
|---|
| 379 | if (INew > IPrev) |
|---|
| 380 | { |
|---|
| 381 | uInsertPos = j; |
|---|
| 382 | IPrev = INew; |
|---|
| 383 | } |
|---|
| 384 | } |
|---|
| 385 | |
|---|
| 386 | // Special case for i=uLengthA, j=uLengthB |
|---|
| 387 | SCORE scoreMax = MPrev[uLengthB-1]; |
|---|
| 388 | int iTraceBack = 0; |
|---|
| 389 | |
|---|
| 390 | SCORE scoreD = DPrev[uLengthB-1] + GapCloseA[uLengthA-1]; |
|---|
| 391 | if (scoreD > scoreMax) |
|---|
| 392 | { |
|---|
| 393 | scoreMax = scoreD; |
|---|
| 394 | iTraceBack = (int) uLengthA - (int) uDeletePos[uLengthB-1]; |
|---|
| 395 | } |
|---|
| 396 | |
|---|
| 397 | SCORE scoreI = IPrev + GapCloseB[uLengthB-1]; |
|---|
| 398 | if (scoreI > scoreMax) |
|---|
| 399 | { |
|---|
| 400 | scoreMax = scoreI; |
|---|
| 401 | iTraceBack = (int) uInsertPos - (int) uLengthB; |
|---|
| 402 | } |
|---|
| 403 | |
|---|
| 404 | TraceBack[uLengthA][uLengthB] = iTraceBack; |
|---|
| 405 | |
|---|
| 406 | TraceBackToPath(TraceBack, uLengthA, uLengthB, Path); |
|---|
| 407 | |
|---|
| 408 | return scoreMax; |
|---|
| 409 | } |
|---|
