9 #if 1 // SINGLE_AFFINE
\r
11 extern bool g_bKeepSimpleDP;
\r
12 extern SCORE *g_DPM;
\r
13 extern SCORE *g_DPD;
\r
14 extern SCORE *g_DPI;
\r
19 static const char *LocalScoreToStr(SCORE s)
\r
21 static char str[16];
\r
24 sprintf(str, "%6.1f", s);
\r
28 static void ListTB(const char *TBM_, const ProfPos *PA, const ProfPos *PB,
\r
29 unsigned uPrefixCountA, unsigned uPrefixCountB)
\r
32 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
35 if (uPrefixLengthB > 0)
\r
36 c = ConsensusChar(PB[uPrefixLengthB - 1]);
\r
37 Log(" %4u:%c", uPrefixLengthB, c);
\r
40 for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
43 if (uPrefixLengthA > 0)
\r
44 c = ConsensusChar(PA[uPrefixLengthA - 1]);
\r
45 Log("%4u:%c ", uPrefixLengthA, c);
\r
46 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
47 Log(" %6c", TBM(uPrefixLengthA, uPrefixLengthB));
\r
52 static void ListDP(const SCORE *DPM_, const ProfPos *PA, const ProfPos *PB,
\r
53 unsigned uPrefixCountA, unsigned uPrefixCountB)
\r
56 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
59 if (uPrefixLengthB > 0)
\r
60 c = ConsensusChar(PB[uPrefixLengthB - 1]);
\r
61 Log(" %4u:%c", uPrefixLengthB, c);
\r
64 for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
67 if (uPrefixLengthA > 0)
\r
68 c = ConsensusChar(PA[uPrefixLengthA - 1]);
\r
69 Log("%4u:%c ", uPrefixLengthA, c);
\r
70 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
71 Log(" %s", LocalScoreToStr(DPM(uPrefixLengthA, uPrefixLengthB)));
\r
76 SCORE GlobalAlignSimple(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,
\r
77 unsigned uLengthB, PWPath &Path)
\r
79 assert(uLengthB > 0 && uLengthA > 0);
\r
81 SetTermGaps(PA, uLengthA);
\r
82 SetTermGaps(PB, uLengthB);
\r
84 const unsigned uPrefixCountA = uLengthA + 1;
\r
85 const unsigned uPrefixCountB = uLengthB + 1;
\r
87 // Allocate DP matrices
\r
88 const size_t LM = uPrefixCountA*uPrefixCountB;
\r
89 SCORE *DPL_ = new SCORE[LM];
\r
90 SCORE *DPM_ = new SCORE[LM];
\r
91 SCORE *DPD_ = new SCORE[LM];
\r
92 SCORE *DPI_ = new SCORE[LM];
\r
94 char *TBM_ = new char[LM];
\r
95 char *TBD_ = new char[LM];
\r
96 char *TBI_ = new char[LM];
\r
98 memset(TBM_, '?', LM);
\r
99 memset(TBD_, '?', LM);
\r
100 memset(TBI_, '?', LM);
\r
103 DPD(0, 0) = MINUS_INFINITY;
\r
104 DPI(0, 0) = MINUS_INFINITY;
\r
106 DPM(1, 0) = MINUS_INFINITY;
\r
107 DPD(1, 0) = PA[0].m_scoreGapOpen;
\r
109 DPI(1, 0) = MINUS_INFINITY;
\r
111 DPM(0, 1) = MINUS_INFINITY;
\r
112 DPD(0, 1) = MINUS_INFINITY;
\r
113 DPI(0, 1) = PB[0].m_scoreGapOpen;
\r
116 // Empty prefix of B is special case
\r
117 for (unsigned uPrefixLengthA = 2; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
119 // M=LetterA+LetterB, impossible with empty prefix
\r
120 DPM(uPrefixLengthA, 0) = MINUS_INFINITY;
\r
123 DPD(uPrefixLengthA, 0) = DPD(uPrefixLengthA - 1, 0) + g_scoreGapExtend;
\r
124 TBD(uPrefixLengthA, 0) = 'D';
\r
126 // I=GapA+LetterB, impossible with empty prefix
\r
127 DPI(uPrefixLengthA, 0) = MINUS_INFINITY;
\r
130 // Empty prefix of A is special case
\r
131 for (unsigned uPrefixLengthB = 2; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
133 // M=LetterA+LetterB, impossible with empty prefix
\r
134 DPM(0, uPrefixLengthB) = MINUS_INFINITY;
\r
136 // D=LetterA+GapB, impossible with empty prefix
\r
137 DPD(0, uPrefixLengthB) = MINUS_INFINITY;
\r
140 DPI(0, uPrefixLengthB) = DPI(0, uPrefixLengthB - 1) + g_scoreGapExtend;
\r
141 TBI(0, uPrefixLengthB) = 'I';
\r
144 // Special case to agree with NWFast, no D-I transitions so...
\r
145 DPD(uLengthA, 0) = MINUS_INFINITY;
\r
146 // DPI(0, uLengthB) = MINUS_INFINITY;
\r
151 SCORE scoreGapCloseB = MINUS_INFINITY;
\r
152 for (unsigned uPrefixLengthB = 1; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
154 const ProfPos &PPB = PB[uPrefixLengthB - 1];
\r
156 SCORE scoreGapCloseA = MINUS_INFINITY;
\r
157 for (unsigned uPrefixLengthA = 1; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
159 const ProfPos &PPA = PA[uPrefixLengthA - 1];
\r
162 // Match M=LetterA+LetterB
\r
163 SCORE scoreLL = ScoreProfPos2(PPA, PPB);
\r
164 DPL(uPrefixLengthA, uPrefixLengthB) = scoreLL;
\r
166 SCORE scoreMM = DPM(uPrefixLengthA-1, uPrefixLengthB-1);
\r
167 SCORE scoreDM = DPD(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseA;
\r
168 SCORE scoreIM = DPI(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseB;
\r
171 if (scoreMM >= scoreDM && scoreMM >= scoreIM)
\r
173 scoreBest = scoreMM;
\r
174 TBM(uPrefixLengthA, uPrefixLengthB) = 'M';
\r
176 else if (scoreDM >= scoreMM && scoreDM >= scoreIM)
\r
178 scoreBest = scoreDM;
\r
179 TBM(uPrefixLengthA, uPrefixLengthB) = 'D';
\r
183 assert(scoreIM >= scoreMM && scoreIM >= scoreDM);
\r
184 scoreBest = scoreIM;
\r
185 TBM(uPrefixLengthA, uPrefixLengthB) = 'I';
\r
187 DPM(uPrefixLengthA, uPrefixLengthB) = scoreBest + scoreLL;
\r
191 // Delete D=LetterA+GapB
\r
192 SCORE scoreMD = DPM(uPrefixLengthA-1, uPrefixLengthB) +
\r
193 PA[uPrefixLengthA-1].m_scoreGapOpen;
\r
194 SCORE scoreDD = DPD(uPrefixLengthA-1, uPrefixLengthB) + g_scoreGapExtend;
\r
197 if (scoreMD >= scoreDD)
\r
199 scoreBest = scoreMD;
\r
200 TBD(uPrefixLengthA, uPrefixLengthB) = 'M';
\r
204 assert(scoreDD >= scoreMD);
\r
205 scoreBest = scoreDD;
\r
206 TBD(uPrefixLengthA, uPrefixLengthB) = 'D';
\r
208 DPD(uPrefixLengthA, uPrefixLengthB) = scoreBest;
\r
211 // Insert I=GapA+LetterB
\r
213 SCORE scoreMI = DPM(uPrefixLengthA, uPrefixLengthB-1) +
\r
214 PB[uPrefixLengthB - 1].m_scoreGapOpen;
\r
215 SCORE scoreII = DPI(uPrefixLengthA, uPrefixLengthB-1) + g_scoreGapExtend;
\r
218 if (scoreMI >= scoreII)
\r
220 scoreBest = scoreMI;
\r
221 TBI(uPrefixLengthA, uPrefixLengthB) = 'M';
\r
225 assert(scoreII > scoreMI);
\r
226 scoreBest = scoreII;
\r
227 TBI(uPrefixLengthA, uPrefixLengthB) = 'I';
\r
229 DPI(uPrefixLengthA, uPrefixLengthB) = scoreBest;
\r
232 scoreGapCloseA = PPA.m_scoreGapClose;
\r
234 scoreGapCloseB = PPB.m_scoreGapClose;
\r
239 Log("Simple DPL:\n");
\r
240 ListDP(DPL_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
242 Log("Simple DPM:\n");
\r
243 ListDP(DPM_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
245 Log("Simple DPD:\n");
\r
246 ListDP(DPD_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
248 Log("Simple DPI:\n");
\r
249 ListDP(DPI_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
251 Log("Simple TBM:\n");
\r
252 ListTB(TBM_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
254 Log("Simple TBD:\n");
\r
255 ListTB(TBD_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
257 Log("Simple TBI:\n");
\r
258 ListTB(TBI_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
266 SCORE M = DPM(uLengthA, uLengthB);
\r
267 SCORE D = DPD(uLengthA, uLengthB) + PA[uLengthA-1].m_scoreGapClose;
\r
268 SCORE I = DPI(uLengthA, uLengthB) + PB[uLengthB-1].m_scoreGapClose;
\r
269 char cEdgeType = '?';
\r
271 SCORE BestScore = MINUS_INFINITY;
\r
272 if (M >= D && M >= I)
\r
277 else if (D >= M && D >= I)
\r
284 assert(I >= M && I >= D);
\r
290 Log("Simple: MAB=%.4g DAB=%.4g IAB=%.4g best=%c\n", M, D, I, cEdgeType);
\r
293 unsigned PLA = uLengthA;
\r
294 unsigned PLB = uLengthB;
\r
298 Edge.cType = cEdgeType;
\r
299 Edge.uPrefixLengthA = PLA;
\r
300 Edge.uPrefixLengthB = PLB;
\r
302 Log("Prepend %c%d.%d\n", Edge.cType, PLA, PLB);
\r
304 Path.PrependEdge(Edge);
\r
311 cEdgeType = TBM(PLA, PLB);
\r
318 cEdgeType = TBD(PLA, PLB);
\r
324 cEdgeType = TBI(PLA, PLB);
\r
329 Quit("Invalid edge %c", cEdgeType);
\r
331 if (0 == PLA && 0 == PLB)
\r
336 // SCORE Score = TraceBack(PA, uLengthA, PB, uLengthB, DPM_, DPD_, DPI_, Path);
\r
339 SCORE scorePath = FastScorePath2(PA, uLengthA, PB, uLengthB, Path);
\r
341 Log("Score = %s Path = %s\n", LocalScoreToStr(BestScore), LocalScoreToStr(scorePath));
\r
344 if (g_bKeepSimpleDP)
\r
368 #endif // SINLGLE_AFFINE
\r