9 S[i] = Score of col i: 0=no SNP, +1 = Y, -3 = N or A.
\r
11 V[k] = Best segment score from j, j+1 .. k for all possible j
\r
12 max(j) Sum i=j..k S[i]
\r
15 V[k] = S[k] + max (V[k-1], 0)
\r
18 void AlignChimeGlobal3(const string &Q3, const string &A3, const string &B3,
\r
19 const string &QLabel, const string &ALabel, const string &BLabel,
\r
22 void Make3Way(const SeqData &SDQ, const SeqData &SDA, const SeqData &SDB,
\r
23 const string &PathQA, const string &PathQB,
\r
24 string &Q3, string &A3, string &B3);
\r
26 double GetScore2(double Y, double N, double A);
\r
28 void AlignChimeLocal3(const string &Q3, const string &A3, const string &B3,
\r
29 const string &QLabel, const string &ALabel, const string &BLabel,
\r
34 const byte *Q3Seq = (const byte *) Q3.c_str();
\r
35 const byte *A3Seq = (const byte *) A3.c_str();
\r
36 const byte *B3Seq = (const byte *) B3.c_str();
\r
38 const unsigned ColCount = SIZE(Q3);
\r
39 asserta(SIZE(A3) == ColCount && SIZE(B3) == ColCount);
\r
41 vector<float> ColScoresA(ColCount, 0.0f);
\r
42 vector<float> ColScoresB(ColCount, 0.0f);
\r
44 float ScoreN = -(float) opt_xn;
\r
46 for (unsigned Col = 0; Col < ColCount; ++Col)
\r
48 char q = Q3Seq[Col];
\r
49 char a = A3Seq[Col];
\r
50 char b = B3Seq[Col];
\r
55 if (q == a && q == b && a == b)
\r
58 if (isgap(q) || isgap(a) || isgap(b))
\r
61 if (Col > 0 && (isgap(Q3Seq[Col-1]) || isgap(A3Seq[Col-1]) || isgap(B3Seq[Col-1])))
\r
64 if (Col + 1 < ColCount && (isgap(Q3Seq[Col+1]) || isgap(A3Seq[Col+1]) || isgap(B3Seq[Col+1])))
\r
67 if (q == a && q != b)
\r
68 ColScoresA[Col] = 1;
\r
70 ColScoresA[Col] = ScoreN;
\r
72 if (q == b && q != a)
\r
73 ColScoresB[Col] = 1;
\r
75 ColScoresB[Col] = ScoreN;
\r
78 vector<float> LVA(ColCount, 0.0f);
\r
79 vector<float> LVB(ColCount, 0.0f);
\r
81 LVA[0] = ColScoresA[0];
\r
82 LVB[0] = ColScoresB[0];
\r
83 for (unsigned Col = 1; Col < ColCount; ++Col)
\r
85 LVA[Col] = max(LVA[Col-1], 0.0f) + ColScoresA[Col];
\r
86 LVB[Col] = max(LVB[Col-1], 0.0f) + ColScoresB[Col];
\r
89 vector<float> RVA(ColCount, 0.0f);
\r
90 vector<float> RVB(ColCount, 0.0f);
\r
92 RVA[ColCount-1] = ColScoresA[ColCount-1];
\r
93 RVB[ColCount-1] = ColScoresB[ColCount-1];
\r
94 for (int Col = ColCount-2; Col >= 0; --Col)
\r
96 RVA[Col] = max(RVA[Col+1], 0.0f) + ColScoresA[Col];
\r
97 RVB[Col] = max(RVB[Col+1], 0.0f) + ColScoresB[Col];
\r
100 bool FirstA = true;
\r
101 float MaxSum = 0.0;
\r
102 unsigned ColX = UINT_MAX;
\r
103 for (unsigned Col = 1; Col < ColCount-1; ++Col)
\r
105 float Sum = LVA[Col] + RVB[Col+1];
\r
114 for (unsigned Col = 1; Col < ColCount-1; ++Col)
\r
116 float Sum = LVB[Col] + RVA[Col+1];
\r
124 if (ColX == UINT_MAX)
\r
127 unsigned ColLo = UINT_MAX;
\r
128 unsigned ColHi = UINT_MAX;
\r
132 for (int Col = ColX; Col >= 0; --Col)
\r
134 Sum += ColScoresA[Col];
\r
135 if (Sum >= LVA[ColX])
\r
141 asserta(Sum >= LVA[ColX]);
\r
143 for (unsigned Col = ColX+1; Col < ColCount; ++Col)
\r
145 Sum += ColScoresB[Col];
\r
146 if (Sum >= RVB[ColX])
\r
152 asserta(Sum >= RVB[ColX]);
\r
157 for (int Col = ColX; Col >= 0; --Col)
\r
159 Sum += ColScoresB[Col];
\r
160 if (Sum >= LVB[ColX])
\r
166 asserta(Sum >= LVB[ColX]);
\r
168 for (unsigned Col = ColX+1; Col < ColCount; ++Col)
\r
170 Sum += ColScoresA[Col];
\r
171 if (Sum >= RVA[ColX])
\r
177 asserta(Sum >= RVA[ColX]);
\r
180 unsigned ColXHi = ColX;
\r
181 for (unsigned Col = ColX + 1; Col < ColCount; ++Col)
\r
183 char q = Q3Seq[Col];
\r
184 char a = A3Seq[Col];
\r
185 char b = B3Seq[Col];
\r
187 if (q == a && q == b && !isgap(q))
\r
193 unsigned ColXLo = ColX;
\r
194 for (int Col = (int) ColX - 1; Col >= 0; --Col)
\r
196 char q = Q3Seq[Col];
\r
197 char a = A3Seq[Col];
\r
198 char b = B3Seq[Col];
\r
200 if (q == a && q == b && !isgap(q))
\r
212 for (unsigned Col = 0; Col < ColCount; ++Col)
\r
214 char q = Q3Seq[Col];
\r
215 char a = A3Seq[Col];
\r
216 char b = B3Seq[Col];
\r
218 if (!isgap(q) && !isgap(a))
\r
225 if (!isgap(q) && !isgap(b))
\r
232 if (!isgap(a) && !isgap(b))
\r
240 Hit.PctIdQA = Pct(IdQA, NQA);
\r
241 Hit.PctIdQB = Pct(IdQB, NQB);
\r
242 Hit.PctIdAB = Pct(IdAB, NAB);
\r
244 unsigned LIdQA = 0;
\r
245 unsigned LIdQB = 0;
\r
246 for (unsigned Col = ColLo; Col < ColXLo; ++Col)
\r
248 char q = Q3Seq[Col];
\r
249 char a = A3Seq[Col];
\r
250 char b = B3Seq[Col];
\r
252 if (!isgap(q) && !isgap(a))
\r
258 if (!isgap(q) && !isgap(b))
\r
265 unsigned RIdQA = 0;
\r
266 unsigned RIdQB = 0;
\r
267 for (unsigned Col = ColXHi+1; Col <= ColHi; ++Col)
\r
269 char q = Q3Seq[Col];
\r
270 char a = A3Seq[Col];
\r
271 char b = B3Seq[Col];
\r
273 if (!isgap(q) && !isgap(a))
\r
279 if (!isgap(q) && !isgap(b))
\r
286 unsigned IdDiffL = max(LIdQA, LIdQB) - min(LIdQA, LIdQB);
\r
287 unsigned IdDiffR = max(RIdQA, RIdQB) - min(RIdQA, RIdQB);
\r
288 unsigned MinIdDiff = min(IdDiffL, IdDiffR);
\r
289 unsigned ColRange = ColHi - ColLo + 1;
\r
290 if (opt_queryfract > 0.0f && float(ColRange)/float(QL) < opt_queryfract)
\r
293 // double Div = Pct(MinIdDiff, QSD.L);
\r
297 Log(" Col A Q B ScoreA ScoreB LVA LVB RVA RVB\n");
\r
298 Log("----- - - - ------- ------- ------- ------- ------- -------\n");
\r
299 for (unsigned Col = 0; Col < ColCount; ++Col)
\r
301 if (ColScoresA[Col] == 0.0 && ColScoresB[Col] == 0.0)
\r
304 char q = Q3Seq[Col];
\r
305 char a = A3Seq[Col];
\r
306 char b = B3Seq[Col];
\r
307 Log("%5u %c %c %c", Col, a, q, b);
\r
309 if (ColScoresA[Col] == 0.0)
\r
312 Log(" %7.1f", ColScoresA[Col]);
\r
314 if (ColScoresB[Col] == 0.0)
\r
317 Log(" %7.1f", ColScoresB[Col]);
\r
319 Log(" %7.1f %7.1f %7.1f %7.1f", LVA[Col], LVB[Col], RVA[Col], RVB[Col]);
\r
324 Log("MaxSum %.1f, ColLo %u, ColXLo %u, ColX %u, ColXHi %u, ColHi %u, AF %c\n",
\r
325 MaxSum, ColLo, ColXLo, ColX, ColXHi, ColHi, tof(FirstA));
\r
326 Log(" LIdQA %u, LIdQB %u, RIdQA %u, RIdQB %u\n", LIdQA, LIdQB, RIdQA, RIdQB);
\r
333 for (unsigned Col = ColLo; Col <= ColHi; ++Col)
\r
344 AlignChimeGlobal3(Q3L, A3L, B3L, QLabel, ALabel, BLabel, Hit);
\r
354 for (unsigned Col = ColLo; Col <= ColHi; ++Col)
\r
356 char q = Q3Seq[Col];
\r
357 char a = A3Seq[Col];
\r
358 char b = B3Seq[Col];
\r
359 if (q == a && q == b && a == b)
\r
361 if (isgap(q) || isgap(a) || isgap(b))
\r
363 if (Col > 0 && (isgap(Q3Seq[Col-1]) || isgap(A3Seq[Col-1]) || isgap(B3Seq[Col-1])))
\r
365 if (Col + 1 < ColCount && (isgap(Q3Seq[Col+1]) || isgap(A3Seq[Col+1]) || isgap(B3Seq[Col+1])))
\r
373 if (q == a && q != b)
\r
375 else if (q == b && q != a)
\r
380 else if (Col > ColXHi)
\r
382 if (q == b && q != a)
\r
384 else if (q == a && q != b)
\r
391 double ScoreL = GetScore2(Hit.CS_LY, Hit.CS_LN, Hit.CS_LA);
\r
392 double ScoreR = GetScore2(Hit.CS_RY, Hit.CS_RN, Hit.CS_RA);
\r
393 Hit.Score = ScoreL*ScoreR;
\r
400 // Hit.PathQA = PathQA;
\r
401 // Hit.PathQB = PathQB;
\r
409 //Hit.ColLo = ColLo;
\r
410 //Hit.ColXLo = ColXLo;
\r
411 //Hit.ColXHi = ColXHi;
\r
412 //Hit.ColHi = ColHi;
\r