source: trunk/GDE/MUSCLE/src/glbalignns.cpp

Last change on this file was 10390, checked in by aboeckma, 12 years ago

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