8 #define EQ(a, b) (fabs(a-b) < 0.1)
\r
10 SCORE TraceBack(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,
\r
11 unsigned uLengthB, const SCORE *DPM_, const SCORE *DPD_, const SCORE *DPI_,
\r
16 Log("TraceBack LengthA=%u LengthB=%u\n", uLengthA, uLengthB);
\r
18 assert(uLengthB > 0 && uLengthA > 0);
\r
20 const unsigned uPrefixCountA = uLengthA + 1;
\r
21 const unsigned uPrefixCountB = uLengthB + 1;
\r
25 unsigned uPrefixLengthA = uLengthA;
\r
26 unsigned uPrefixLengthB = uLengthB;
\r
28 const SCORE scoreM = DPM(uPrefixLengthA, uPrefixLengthB);
\r
29 SCORE scoreD = DPD(uPrefixLengthA, uPrefixLengthB);
\r
30 SCORE scoreI = DPI(uPrefixLengthA, uPrefixLengthB);
\r
32 const ProfPos &LastPPA = PA[uLengthA - 1];
\r
33 const ProfPos &LastPPB = PB[uLengthB - 1];
\r
35 scoreD += LastPPA.m_scoreGapClose;
\r
36 scoreI += LastPPB.m_scoreGapClose;
\r
38 char cEdgeType = cInsane;
\r
40 if (scoreM >= scoreD && scoreM >= scoreI)
\r
45 else if (scoreD >= scoreM && scoreD >= scoreI)
\r
52 assert(scoreI >= scoreM && scoreI >= scoreD);
\r
59 if ('S' == cEdgeType)
\r
63 Edge.cType = cEdgeType;
\r
64 Edge.uPrefixLengthA = uPrefixLengthA;
\r
65 Edge.uPrefixLengthB = uPrefixLengthB;
\r
66 Path.PrependEdge(Edge);
\r
69 unsigned uPrevPrefixLengthA = uPrefixLengthA;
\r
70 unsigned uPrevPrefixLengthB = uPrefixLengthB;
\r
76 assert(uPrefixLengthA > 0);
\r
77 assert(uPrefixLengthB > 0);
\r
78 const ProfPos &PPA = PA[uPrefixLengthA - 1];
\r
79 const ProfPos &PPB = PB[uPrefixLengthB - 1];
\r
81 const SCORE Score = DPM(uPrefixLengthA, uPrefixLengthB);
\r
82 const SCORE scoreMatch = ScoreProfPos2(PPA, PPB);
\r
85 if (1 == uPrefixLengthA && 1 == uPrefixLengthB)
\r
86 scoreSM = scoreMatch;
\r
88 scoreSM = MINUS_INFINITY;
\r
90 SCORE scoreMM = MINUS_INFINITY;
\r
91 SCORE scoreDM = MINUS_INFINITY;
\r
92 SCORE scoreIM = MINUS_INFINITY;
\r
93 if (uPrefixLengthA > 1 && uPrefixLengthB > 1)
\r
94 scoreMM = DPM(uPrefixLengthA-1, uPrefixLengthB-1) + scoreMatch;
\r
95 if (uPrefixLengthA > 1)
\r
97 SCORE scoreTransDM = PA[uPrefixLengthA-2].m_scoreGapClose;
\r
98 scoreDM = DPD(uPrefixLengthA-1, uPrefixLengthB-1) + scoreTransDM + scoreMatch;
\r
100 if (uPrefixLengthB > 1)
\r
102 SCORE scoreTransIM = PB[uPrefixLengthB-2].m_scoreGapClose;
\r
103 scoreIM = DPI(uPrefixLengthA-1, uPrefixLengthB-1) + scoreTransIM + scoreMatch;
\r
106 if (EQ(scoreMM, Score))
\r
107 cPrevEdgeType = 'M';
\r
108 else if (EQ(scoreDM, Score))
\r
109 cPrevEdgeType = 'D';
\r
110 else if (EQ(scoreIM, Score))
\r
111 cPrevEdgeType = 'I';
\r
112 else if (EQ(scoreSM, Score))
\r
113 cPrevEdgeType = 'S';
\r
115 Quit("TraceBack: failed to match M score=%g M=%g D=%g I=%g S=%g",
\r
116 Score, scoreMM, scoreDM, scoreIM, scoreSM);
\r
118 --uPrevPrefixLengthA;
\r
119 --uPrevPrefixLengthB;
\r
125 assert(uPrefixLengthA > 0);
\r
126 const SCORE Score = DPD(uPrefixLengthA, uPrefixLengthB);
\r
128 SCORE scoreMD = MINUS_INFINITY;
\r
129 SCORE scoreDD = MINUS_INFINITY;
\r
130 SCORE scoreSD = MINUS_INFINITY;
\r
131 if (uPrefixLengthB == 0)
\r
133 if (uPrefixLengthA == 1)
\r
134 scoreSD = PA[0].m_scoreGapOpen;
\r
136 scoreSD = DPD(uPrefixLengthA - 1, 0);
\r
138 if (uPrefixLengthA > 1)
\r
140 const ProfPos &PPA = PA[uPrefixLengthA - 1];
\r
141 SCORE scoreTransMD = PPA.m_scoreGapOpen;
\r
142 scoreMD = DPM(uPrefixLengthA-1, uPrefixLengthB) + scoreTransMD;
\r
143 scoreDD = DPD(uPrefixLengthA-1, uPrefixLengthB);
\r
146 if (EQ(Score, scoreMD))
\r
147 cPrevEdgeType = 'M';
\r
148 else if (EQ(Score, scoreDD))
\r
149 cPrevEdgeType = 'D';
\r
150 else if (EQ(Score, scoreSD))
\r
151 cPrevEdgeType = 'S';
\r
153 Quit("TraceBack: failed to match D");
\r
155 --uPrevPrefixLengthA;
\r
161 assert(uPrefixLengthB > 0);
\r
162 const SCORE Score = DPI(uPrefixLengthA, uPrefixLengthB);
\r
164 SCORE scoreMI = MINUS_INFINITY;
\r
165 SCORE scoreII = MINUS_INFINITY;
\r
166 SCORE scoreSI = MINUS_INFINITY;
\r
167 if (uPrefixLengthA == 0)
\r
169 if (uPrefixLengthB == 1)
\r
170 scoreSI = PB[0].m_scoreGapOpen;
\r
172 scoreSI = DPI(0, uPrefixLengthB - 1);
\r
174 if (uPrefixLengthB > 1)
\r
176 const ProfPos &PPB = PB[uPrefixLengthB - 1];
\r
177 SCORE scoreTransMI = PPB.m_scoreGapOpen;
\r
178 scoreMI = DPM(uPrefixLengthA, uPrefixLengthB-1) + scoreTransMI;
\r
179 scoreII = DPI(uPrefixLengthA, uPrefixLengthB-1);
\r
182 if (EQ(Score, scoreMI))
\r
183 cPrevEdgeType = 'M';
\r
184 else if (EQ(Score, scoreII))
\r
185 cPrevEdgeType = 'I';
\r
186 else if (EQ(Score, scoreSI))
\r
187 cPrevEdgeType = 'S';
\r
189 Quit("TraceBack: failed to match I");
\r
191 --uPrevPrefixLengthB;
\r
199 Log("Edge %c%c%u.%u", cPrevEdgeType, cEdgeType, uPrefixLengthA, uPrefixLengthB);
\r
202 cEdgeType = cPrevEdgeType;
\r
203 uPrefixLengthA = uPrevPrefixLengthA;
\r
204 uPrefixLengthB = uPrevPrefixLengthB;
\r