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
FreeDPMemSPN()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
AllocDPMem(unsigned uLengthA,unsigned uLengthB)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
GlobalAlignSPN(const ProfPos * PA,unsigned uLengthA,const ProfPos * PB,unsigned uLengthB,PWPath & Path)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 }
410