3 #include <stdio.h> // for sprintf
\r
6 #include "gapscoredimer.h"
\r
10 static SCORE TraceBackDimer( const SCORE *DPM_, const SCORE *DPD_, const SCORE *DPI_,
\r
11 const char *TBM_, const char *TBD_, const char *TBI_,
\r
12 unsigned uLengthA, unsigned uLengthB, PWPath &Path);
\r
14 static const char *LocalScoreToStr(SCORE s)
\r
16 static char str[16];
\r
17 if (MINUS_INFINITY == s)
\r
19 sprintf(str, "%6.3g", s);
\r
24 static void ListDP(const SCORE *DPM_, const ProfPos *PA, const ProfPos *PB,
\r
25 unsigned uPrefixCountA, unsigned uPrefixCountB)
\r
28 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
31 if (uPrefixLengthB > 0)
\r
32 c = ConsensusChar(PB[uPrefixLengthB - 1]);
\r
33 Log(" %4u:%c", uPrefixLengthB, c);
\r
36 for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
39 if (uPrefixLengthA > 0)
\r
40 c = ConsensusChar(PA[uPrefixLengthA - 1]);
\r
41 Log("%4u:%c ", uPrefixLengthA, c);
\r
42 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
43 Log(" %s", LocalScoreToStr(DPM(uPrefixLengthA, uPrefixLengthB)));
\r
48 static void ListTB(const char *TBM_, const ProfPos *PA, const ProfPos *PB,
\r
49 unsigned uPrefixCountA, unsigned uPrefixCountB)
\r
52 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
53 Log("%2d", uPrefixLengthB);
\r
56 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
59 if (uPrefixLengthB > 0)
\r
60 c = ConsensusChar(PB[uPrefixLengthB - 1]);
\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(" %c", TBM(uPrefixLengthA, uPrefixLengthB));
\r
77 static ProfPos PPTerm;
\r
78 static bool InitializePPTerm()
\r
80 PPTerm.m_bAllGaps = false;
\r
88 static bool PPTermInitialized = InitializePPTerm();
\r
90 static SCORE ScoreProfPosDimerLE(const ProfPos &PPA, const ProfPos &PPB)
\r
93 for (unsigned n = 0; n < 20; ++n)
\r
95 const unsigned uLetter = PPA.m_uSortOrder[n];
\r
96 const FCOUNT fcLetter = PPA.m_fcCounts[uLetter];
\r
99 Score += fcLetter*PPB.m_AAScores[uLetter];
\r
103 SCORE logScore = logf(Score);
\r
104 return (SCORE) (logScore*(PPA.m_fOcc * PPB.m_fOcc));
\r
107 static SCORE ScoreProfPosDimerPSP(const ProfPos &PPA, const ProfPos &PPB)
\r
110 for (unsigned n = 0; n < 20; ++n)
\r
112 const unsigned uLetter = PPA.m_uSortOrder[n];
\r
113 const FCOUNT fcLetter = PPA.m_fcCounts[uLetter];
\r
116 Score += fcLetter*PPB.m_AAScores[uLetter];
\r
121 static SCORE ScoreProfPosDimer(const ProfPos &PPA, const ProfPos &PPB)
\r
126 return ScoreProfPosDimerLE(PPA, PPB);
\r
130 return ScoreProfPosDimerPSP(PPA, PPB);
\r
132 Quit("Invalid g_PPScore");
\r
136 // Global alignment dynamic programming
\r
137 // This variant optimizes the profile-profile SP score under the
\r
138 // dimer approximation.
\r
139 SCORE GlobalAlignDimer(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,
\r
140 unsigned uLengthB, PWPath &Path)
\r
142 assert(uLengthB > 0 && uLengthA > 0);
\r
144 const unsigned uPrefixCountA = uLengthA + 1;
\r
145 const unsigned uPrefixCountB = uLengthB + 1;
\r
147 // Allocate DP matrices
\r
148 const size_t LM = uPrefixCountA*uPrefixCountB;
\r
149 SCORE *DPM_ = new SCORE[LM];
\r
150 SCORE *DPD_ = new SCORE[LM];
\r
151 SCORE *DPI_ = new SCORE[LM];
\r
153 char *TBM_ = new char[LM];
\r
154 char *TBD_ = new char[LM];
\r
155 char *TBI_ = new char[LM];
\r
158 DPD(0, 0) = MINUS_INFINITY;
\r
159 DPI(0, 0) = MINUS_INFINITY;
\r
165 DPM(1, 0) = MINUS_INFINITY;
\r
166 DPD(1, 0) = GapScoreMD(PA[0], PPTerm);
\r
167 DPI(1, 0) = MINUS_INFINITY;
\r
173 DPM(0, 1) = MINUS_INFINITY;
\r
174 DPD(0, 1) = MINUS_INFINITY;
\r
175 DPI(0, 1) = GapScoreMI(PPTerm, PB[0]);
\r
181 // Empty prefix of B is special case
\r
182 for (unsigned uPrefixLengthA = 2; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
184 // M=LetterA+LetterB, impossible with empty prefix
\r
185 DPM(uPrefixLengthA, 0) = MINUS_INFINITY;
\r
186 TBM(uPrefixLengthA, 0) = '?';
\r
189 DPD(uPrefixLengthA, 0) = DPD(uPrefixLengthA - 1, 0) +
\r
190 GapScoreDD(PA[uPrefixLengthA - 1], PPTerm);
\r
191 TBD(uPrefixLengthA, 0) = 'D';
\r
193 // I=GapA+LetterB, impossible with empty prefix
\r
194 DPI(uPrefixLengthA, 0) = MINUS_INFINITY;
\r
195 TBI(uPrefixLengthA, 0) = '?';
\r
198 // Empty prefix of A is special case
\r
199 for (unsigned uPrefixLengthB = 2; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
201 // M=LetterA+LetterB, impossible with empty prefix
\r
202 DPM(0, uPrefixLengthB) = MINUS_INFINITY;
\r
203 TBM(0, uPrefixLengthB) = '?';
\r
205 // D=LetterA+GapB, impossible with empty prefix
\r
206 DPD(0, uPrefixLengthB) = MINUS_INFINITY;
\r
207 TBD(0, uPrefixLengthB) = '?';
\r
210 DPI(0, uPrefixLengthB) = DPI(0, uPrefixLengthB - 1) +
\r
211 GapScoreII(PPTerm, PB[uPrefixLengthB - 1]);
\r
212 TBI(0, uPrefixLengthB) = 'I';
\r
218 for (unsigned uPrefixLengthB = 1; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
220 const ProfPos &PPB = PB[uPrefixLengthB - 1];
\r
221 for (unsigned uPrefixLengthA = 1; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
223 const ProfPos &PPA = PA[uPrefixLengthA - 1];
\r
225 // Match M=LetterA+LetterB
\r
226 SCORE scoreLL = ScoreProfPosDimer(PPA, PPB);
\r
228 SCORE scoreMM = DPM(uPrefixLengthA-1, uPrefixLengthB-1) + GapScoreMM(PPA, PPB);
\r
229 SCORE scoreDM = DPD(uPrefixLengthA-1, uPrefixLengthB-1) + GapScoreDM(PPA, PPB);
\r
230 SCORE scoreIM = DPI(uPrefixLengthA-1, uPrefixLengthB-1) + GapScoreIM(PPA, PPB);
\r
232 SCORE scoreBest = scoreMM;
\r
234 if (scoreDM > scoreBest)
\r
236 scoreBest = scoreDM;
\r
239 if (scoreIM > scoreBest)
\r
241 scoreBest = scoreIM;
\r
245 DPM(uPrefixLengthA, uPrefixLengthB) = scoreBest + scoreLL;
\r
246 TBM(uPrefixLengthA, uPrefixLengthB) = c;
\r
249 // Delete D=LetterA+GapB
\r
250 SCORE scoreMD = DPM(uPrefixLengthA-1, uPrefixLengthB) + GapScoreMD(PPA, PPB);
\r
251 SCORE scoreDD = DPD(uPrefixLengthA-1, uPrefixLengthB) + GapScoreDD(PPA, PPB);
\r
252 SCORE scoreID = DPI(uPrefixLengthA-1, uPrefixLengthB) + GapScoreID(PPA, PPB);
\r
254 SCORE scoreBest = scoreMD;
\r
256 if (scoreDD > scoreBest)
\r
258 scoreBest = scoreDD;
\r
261 if (scoreID > scoreBest)
\r
263 scoreBest = scoreID;
\r
267 DPD(uPrefixLengthA, uPrefixLengthB) = scoreBest;
\r
268 TBD(uPrefixLengthA, uPrefixLengthB) = c;
\r
271 // Insert I=GapA+LetterB
\r
272 SCORE scoreMI = DPM(uPrefixLengthA, uPrefixLengthB-1) + GapScoreMI(PPA, PPB);
\r
273 SCORE scoreDI = DPD(uPrefixLengthA, uPrefixLengthB-1) + GapScoreDI(PPA, PPB);
\r
274 SCORE scoreII = DPI(uPrefixLengthA, uPrefixLengthB-1) + GapScoreII(PPA, PPB);
\r
276 SCORE scoreBest = scoreMI;
\r
278 if (scoreDI > scoreBest)
\r
280 scoreBest = scoreDI;
\r
283 if (scoreII > scoreBest)
\r
285 scoreBest = scoreII;
\r
289 DPI(uPrefixLengthA, uPrefixLengthB) = scoreBest;
\r
290 TBI(uPrefixLengthA, uPrefixLengthB) = c;
\r
297 ListDP(DPM_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
299 ListDP(DPD_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
301 ListDP(DPI_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
303 ListTB(TBM_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
305 ListTB(TBD_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
307 ListTB(TBI_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
310 SCORE Score = TraceBackDimer(DPM_, DPD_, DPI_, TBM_, TBD_, TBI_,
\r
311 uLengthA, uLengthB, Path);
\r
314 Log("GlobalAlignDimer score = %.3g\n", Score);
\r
328 static SCORE TraceBackDimer( const SCORE *DPM_, const SCORE *DPD_, const SCORE *DPI_,
\r
329 const char *TBM_, const char *TBD_, const char *TBI_,
\r
330 unsigned uLengthA, unsigned uLengthB, PWPath &Path)
\r
332 const unsigned uPrefixCountA = uLengthA + 1;
\r
334 unsigned uPrefixLengthA = uLengthA;
\r
335 unsigned uPrefixLengthB = uLengthB;
\r
338 SCORE scoreMax = DPM(uLengthA, uLengthB);
\r
339 if (DPD(uLengthA, uLengthB) > scoreMax)
\r
341 scoreMax = DPD(uLengthA, uLengthB);
\r
344 if (DPI(uLengthA, uLengthB) > scoreMax)
\r
346 scoreMax = DPI(uLengthA, uLengthB);
\r
352 if (0 == uPrefixLengthA && 0 == uPrefixLengthB)
\r
356 Edge.cType = cEdge;
\r
357 Edge.uPrefixLengthA = uPrefixLengthA;
\r
358 Edge.uPrefixLengthB = uPrefixLengthB;
\r
359 Path.PrependEdge(Edge);
\r
362 Log("PLA=%u PLB=%u Edge=%c\n", uPrefixLengthA, uPrefixLengthB, cEdge);
\r
367 assert(uPrefixLengthA > 0 && uPrefixLengthB > 0);
\r
368 cEdge = TBM(uPrefixLengthA, uPrefixLengthB);
\r
373 assert(uPrefixLengthA > 0);
\r
374 cEdge = TBD(uPrefixLengthA, uPrefixLengthB);
\r
378 assert(uPrefixLengthB > 0);
\r
379 cEdge = TBI(uPrefixLengthA, uPrefixLengthB);
\r
383 Quit("Invalid edge PLA=%u PLB=%u %c", uPrefixLengthA, uPrefixLengthB, cEdge);
\r