6 extern SCOREMATRIX VTML_SP;
\r
8 // #define SUBST(i, j) Subst(seqA, seqB, i, j)
\r
9 #define SUBST(i, j) MxRowA[i][seqB.GetLetter(j)]
\r
11 static SCORE Subst(const Seq &seqA, const Seq &seqB, unsigned i, unsigned j)
\r
13 assert(i < seqA.Length());
\r
14 assert(j < seqB.Length());
\r
16 unsigned uLetterA = seqA.GetLetter(i);
\r
17 unsigned uLetterB = seqB.GetLetter(j);
\r
18 return VTML_SP[uLetterA][uLetterB] + g_scoreCenter;
\r
32 unsigned *uDeletePos;
\r
36 static struct DP_MEMORY DPM;
\r
38 static void AllocDPMem(unsigned uLengthA, unsigned uLengthB)
\r
40 // Max prefix length
\r
41 unsigned uLength = (uLengthA > uLengthB ? uLengthA : uLengthB) + 1;
\r
42 if (uLength < DPM.uLength)
\r
45 // Add 256 to allow for future expansion and
\r
46 // round up to next multiple of 32.
\r
48 uLength += 32 - uLength%32;
\r
50 const unsigned uOldLength = DPM.uLength;
\r
53 for (unsigned i = 0; i < uOldLength; ++i)
\r
54 delete[] DPM.TraceBack[i];
\r
62 delete[] DPM.MxRowA;
\r
63 delete[] DPM.LettersB;
\r
64 delete[] DPM.uDeletePos;
\r
65 delete[] DPM.TraceBack;
\r
68 DPM.uLength = uLength;
\r
70 DPM.MPrev = new SCORE[uLength];
\r
71 DPM.MCurr = new SCORE[uLength];
\r
72 DPM.MWork = new SCORE[uLength];
\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
81 DPM.TraceBack = new int*[uLength];
\r
83 for (unsigned i = 0; i < uLength; ++i)
\r
84 DPM.TraceBack[i] = new int[uLength];
\r
87 static void RowFromSeq(const Seq &s, SCORE *Row[])
\r
89 const unsigned uLength = s.Length();
\r
90 for (unsigned i = 0; i < uLength; ++i)
\r
92 char c = s.GetChar(i);
\r
93 unsigned uLetter = CharToLetter(c);
\r
95 Row[i] = VTML_SP[uLetter];
\r
97 Row[i] = VTML_SP[AX_X];
\r
101 static void LettersFromSeq(const Seq &s, unsigned Letters[])
\r
103 const unsigned uLength = s.Length();
\r
104 for (unsigned i = 0; i < uLength; ++i)
\r
106 char c = s.GetChar(i);
\r
107 unsigned uLetter = CharToLetter(c);
\r
109 Letters[i] = uLetter;
\r
115 SCORE GlobalAlignSS(const Seq &seqA, const Seq &seqB, PWPath &Path)
\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
122 AllocDPMem(uLengthA, uLengthB);
\r
124 SCORE *MPrev = DPM.MPrev;
\r
125 SCORE *MCurr = DPM.MCurr;
\r
126 SCORE *MWork = DPM.MWork;
\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
134 RowFromSeq(seqA, MxRowA);
\r
135 LettersFromSeq(seqB, LettersB);
\r
137 unsigned *uDeletePos = DPM.uDeletePos;
\r
139 int **TraceBack = DPM.TraceBack;
\r
142 for (unsigned i = 0; i < uPrefixCountA; ++i)
\r
143 memset(TraceBack[i], 0, uPrefixCountB*sizeof(int));
\r
146 // Special case for i=0
\r
147 TraceBack[0][0] = 0;
\r
148 MPrev[0] = MxRowA[0][LettersB[0]];
\r
150 // D(0,0) is -infinity (requires I->D).
\r
151 DPrev[0] = MINUS_INFINITY;
\r
153 for (unsigned j = 1; j < uLengthB; ++j)
\r
155 unsigned uLetterB = LettersB[j];
\r
157 // Only way to get M(0, j) looks like this:
\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
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
171 for (unsigned i = 1; i < uLengthA; ++i)
\r
173 SCORE *ptrMCurr_j = MCurr;
\r
174 memset(ptrMCurr_j, 0, uLengthB*sizeof(SCORE));
\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
182 *ptrMCurr_j = RowA[*ptrLettersB];
\r
186 unsigned *ptrDeletePos = uDeletePos;
\r
188 // Special case for j=0
\r
189 // Only way to get M(i, 0) looks like this:
\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
200 int *ptrTraceBack_ij = TraceBack[i];
\r
201 *ptrTraceBack_ij++ = (int) i;
\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
213 SCORE *ptrDCurr = DCurr;
\r
215 assert(ptrDCurr == &(DCurr[0]));
\r
218 // Can't have an insert if no letters from B
\r
219 IPrev_j_1 = MINUS_INFINITY;
\r
221 unsigned uInsertPos;
\r
223 for (unsigned j = 1; j < uLengthB; ++j)
\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
235 SCORE scoreMax = MPrev_j;
\r
237 assert(ptrDPrev == &(DPrev[j-1]));
\r
238 SCORE scoreD = *ptrDPrev++;
\r
239 if (scoreD > scoreMax)
\r
242 assert(ptrDeletePos == &(uDeletePos[j-1]));
\r
243 *ptrTraceBack_ij = (int) i - (int) *ptrDeletePos;
\r
244 assert(*ptrTraceBack_ij > 0);
\r
248 SCORE scoreI = IPrev_j_1;
\r
249 if (scoreI > scoreMax)
\r
252 *ptrTraceBack_ij = (int) uInsertPos - (int) j;
\r
253 assert(*ptrTraceBack_ij < 0);
\r
256 *ptrMCurr_j += scoreMax;
\r
257 assert(ptrMCurr_j == &(MCurr[j]));
\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
267 assert(ptrDeletePos == &uDeletePos[j]);
\r
270 assert(ptrDCurr + 1 == &(DCurr[j]));
\r
276 Rotate(MPrev, MCurr, MWork);
\r
277 Rotate(DPrev, DCurr, DWork);
\r
280 // Special case for i=uLengthA
\r
281 SCORE IPrev = MINUS_INFINITY;
\r
283 unsigned uInsertPos;
\r
285 for (unsigned j = 1; j < uLengthB; ++j)
\r
287 SCORE INew = MPrev[j-1];
\r
295 // Special case for i=uLengthA, j=uLengthB
\r
296 SCORE scoreMax = MPrev[uLengthB-1];
\r
297 int iTraceBack = 0;
\r
299 SCORE scoreD = DPrev[uLengthB-1] - g_scoreGapOpen/2; // term gaps half
\r
300 if (scoreD > scoreMax)
\r
303 iTraceBack = (int) uLengthA - (int) uDeletePos[uLengthB-1];
\r
306 SCORE scoreI = IPrev - g_scoreGapOpen/2;
\r
307 if (scoreI > scoreMax)
\r
310 iTraceBack = (int) uInsertPos - (int) uLengthB;
\r
313 TraceBack[uLengthA][uLengthB] = iTraceBack;
\r
315 TraceBackToPath(TraceBack, uLengthA, uLengthB, Path);
\r