1 #include "myutils.h"
2 #include "chime.h"
3
WriteChimeFileHdr(FILE * f)4 void WriteChimeFileHdr(FILE *f)
5 {
6 if (f == 0)
7 return;
8
9 fprintf(f,
10 "\tQuery" // 1
11 "\tA" // 2
12 "\tB" // 3
13 "\tIdQM" // 4
14 "\tIdQA" // 5
15 "\tIdQB" // 6
16 "\tIdAB" // 7
17 "\tIdQT" // 8
18 "\tLY" // 9
19 "\tLN" // 10
20 "\tLA" // 11
21 "\tRY" // 12
22 "\tRN" // 13
23 "\tRA" // 14
24 "\tDiv" // 15
25 "\tY" // 16
26 "\n"
27 );
28 }
29
WriteChimeHit(FILE * f,const ChimeHit2 & Hit)30 void WriteChimeHit(FILE *f, const ChimeHit2 &Hit)
31 {
32 if (f == 0)
33 return;
34
35 if (Hit.Div <= 0.0)
36 {
37 fprintf(f, "0.0000"); // 0
38
39 fprintf(f,
40 "\t%s", Hit.QLabel.c_str()); // 1
41
42 fprintf(f,
43 "\t*" // 2
44 "\t*" // 3
45 "\t*" // 4
46 "\t*" // 5
47 "\t*" // 6
48 "\t*" // 7
49 "\t*" // 8
50 "\t*" // 9
51 "\t*" // 10
52 "\t*" // 11
53 "\t*" // 12
54 "\t*" // 13
55 "\t*" // 14
56 "\t*" // 15
57 "\tN" // 16
58 "\n"
59 );
60 return;
61 }
62
63 fprintf(f, "%.4f", Hit.Score); // 0
64
65 fputc('\t', f);
66 fputs(Hit.QLabel.c_str(), f); // 1
67
68 fputc('\t', f);
69 fputs(Hit.ALabel.c_str(), f); // 2
70
71 fputc('\t', f);
72 fputs(Hit.BLabel.c_str(), f); // 3
73
74 fprintf(f, "\t%.1f", Hit.PctIdQM); // 4
75 fprintf(f, "\t%.1f", Hit.PctIdQA); // 5
76 fprintf(f, "\t%.1f", Hit.PctIdQB); // 6
77 fprintf(f, "\t%.1f", Hit.PctIdAB); // 7
78 fprintf(f, "\t%.1f", Hit.PctIdQT); // 8
79
80 fprintf(f, "\t%u", Hit.CS_LY); // 9
81 fprintf(f, "\t%u", Hit.CS_LN); // 10
82 fprintf(f, "\t%u", Hit.CS_LA); // 11
83
84 fprintf(f, "\t%u", Hit.CS_RY); // 12
85 fprintf(f, "\t%u", Hit.CS_RN); // 13
86 fprintf(f, "\t%u", Hit.CS_RA); // 14
87
88 fprintf(f, "\t%.2f", Hit.Div); // 15
89
90 fprintf(f, "\t%c", yon(Hit.Accept())); // 16
91 fputc('\n', f);
92 }
93
GetUngappedLength(const byte * Seq,unsigned L)94 unsigned GetUngappedLength(const byte *Seq, unsigned L)
95 {
96 unsigned UL = 0;
97 for (unsigned i = 0; i < L; ++i)
98 if (!isgap(Seq[i]))
99 ++UL;
100 return UL;
101 }
102
WriteChimeHitX(FILE * f,const ChimeHit2 & Hit)103 void WriteChimeHitX(FILE *f, const ChimeHit2 &Hit)
104 {
105 if (f == 0)
106 return;
107
108 if (Hit.Div <= 0.0)
109 return;
110
111 const string &Q3 = Hit.Q3;
112 const string &A3 = Hit.A3;
113 const string &B3 = Hit.B3;
114
115 const byte *Q3Seq = (const byte *) Q3.c_str();
116 const byte *A3Seq = (const byte *) A3.c_str();
117 const byte *B3Seq = (const byte *) B3.c_str();
118
119 // Aligned
120 unsigned ColCount = SIZE(Q3);
121 asserta(SIZE(A3) == ColCount && SIZE(B3) == ColCount);
122
123 unsigned LQ = GetUngappedLength(Q3Seq, ColCount);
124 unsigned LA = GetUngappedLength(A3Seq, ColCount);
125 unsigned LB = GetUngappedLength(B3Seq, ColCount);
126
127 fprintf(f, "\n");
128 fprintf(f, "------------------------------------------------------------------------\n");
129 fprintf(f, "Query (%5u nt) %s\n", LQ, Hit.QLabel.c_str());
130 fprintf(f, "ParentA (%5u nt) %s\n", LA, Hit.ALabel.c_str());
131 fprintf(f, "ParentB (%5u nt) %s\n", LB, Hit.BLabel.c_str());
132
133 // Strip terminal gaps in query
134 unsigned FromCol = UINT_MAX;
135 unsigned ToCol = UINT_MAX;
136 for (unsigned Col = 0; Col < ColCount; ++Col)
137 {
138 if (!isgap(Q3Seq[Col]))
139 {
140 if (FromCol == UINT_MAX)
141 FromCol = Col;
142 ToCol = Col;
143 }
144 }
145
146 unsigned QPos = 0;
147 unsigned APos = 0;
148 unsigned BPos = 0;
149 for (unsigned Col = 0; Col < FromCol; ++Col)
150 {
151 if (!isgap(A3Seq[Col]))
152 ++APos;
153 if (!isgap(B3Seq[Col]))
154 ++BPos;
155 }
156
157 unsigned Range = ToCol - FromCol + 1;
158 unsigned RowCount = (Range + 79)/80;
159 unsigned RowFromCol = FromCol;
160 for (unsigned RowIndex = 0; RowIndex < RowCount; ++RowIndex)
161 {
162 fprintf(f, "\n");
163 unsigned RowToCol = RowFromCol + 79;
164 if (RowToCol > ToCol)
165 RowToCol = ToCol;
166
167 // A row
168 fprintf(f, "A %5u ", APos + 1);
169 for (unsigned Col = RowFromCol; Col <= RowToCol; ++Col)
170 {
171 char q = Q3Seq[Col];
172 char a = A3Seq[Col];
173 if (a != q)
174 a = tolower(a);
175 fprintf(f, "%c", a);
176 if (!isgap(a))
177 ++APos;
178 }
179 fprintf(f, " %u\n", APos);
180
181 // Q row
182 fprintf(f, "Q %5u ", QPos + 1);
183 for (unsigned Col = RowFromCol; Col <= RowToCol; ++Col)
184 {
185 char q = Q3Seq[Col];
186 fprintf(f, "%c", q);
187 if (!isgap(q))
188 ++QPos;
189 }
190 fprintf(f, " %u\n", QPos);
191
192 // B row
193 fprintf(f, "B %5u ", BPos + 1);
194 for (unsigned Col = RowFromCol; Col <= RowToCol; ++Col)
195 {
196 char q = Q3Seq[Col];
197 char b = B3Seq[Col];
198 if (b != q)
199 b = tolower(b);
200 fprintf(f, "%c", b);
201 if (!isgap(b))
202 ++BPos;
203 }
204 fprintf(f, " %u\n", BPos);
205
206 // Diffs
207 fprintf(f, "Diffs ");
208 for (unsigned Col = RowFromCol; Col <= RowToCol; ++Col)
209 {
210 char q = Q3Seq[Col];
211 char a = A3Seq[Col];
212 char b = B3Seq[Col];
213
214 char c = ' ';
215 if (isgap(q) || isgap(a) || isgap(b))
216 c = ' ';
217 else if (Col < Hit.ColXLo)
218 {
219 if (q == a && q == b)
220 c = ' ';
221 else if (q == a && q != b)
222 c = 'A';
223 else if (q == b && q != a)
224 c = 'b';
225 else if (a == b && q != a)
226 c = 'N';
227 else
228 c = '?';
229 }
230 else if (Col > Hit.ColXHi)
231 {
232 if (q == a && q == b)
233 c = ' ';
234 else if (q == b && q != a)
235 c = 'B';
236 else if (q == a && q != b)
237 c = 'a';
238 else if (a == b && q != a)
239 c = 'N';
240 else
241 c = '?';
242 }
243
244 fprintf(f, "%c", c);
245 }
246 fprintf(f, "\n");
247
248 // SNPs
249 fprintf(f, "Votes ");
250 for (unsigned Col = RowFromCol; Col <= RowToCol; ++Col)
251 {
252 char q = Q3Seq[Col];
253 char a = A3Seq[Col];
254 char b = B3Seq[Col];
255
256 bool PrevGap = Col > 0 && (isgap(Q3Seq[Col-1]) || isgap(A3Seq[Col-1]) || isgap(B3Seq[Col-1]));
257 bool NextGap = Col+1 < ColCount && (isgap(Q3Seq[Col+1]) || isgap(A3Seq[Col+1]) || isgap(B3Seq[Col+1]));
258
259 char c = ' ';
260 if (isgap(q) || isgap(a) || isgap(b) || PrevGap || NextGap)
261 c = ' ';
262 else if (Col < Hit.ColXLo)
263 {
264 if (q == a && q == b)
265 c = ' ';
266 else if (q == a && q != b)
267 c = '+';
268 else if (q == b && q != a)
269 c = '!';
270 else
271 c = '0';
272 }
273 else if (Col > Hit.ColXHi)
274 {
275 if (q == a && q == b)
276 c = ' ';
277 else if (q == b && q != a)
278 c = '+';
279 else if (q == a && q != b)
280 c = '!';
281 else
282 c = '0';
283 }
284
285 fprintf(f, "%c", c);
286 }
287 fprintf(f, "\n");
288
289 // LR row
290 fprintf(f, "Model ");
291 for (unsigned Col = RowFromCol; Col <= RowToCol; ++Col)
292 {
293 if (Col < Hit.ColXLo)
294 fprintf(f, "A");
295 else if (Col >= Hit.ColXLo && Col <= Hit.ColXHi)
296 fprintf(f, "x");
297 else
298 fprintf(f, "B");
299 }
300
301 fprintf(f, "\n");
302
303 RowFromCol += 80;
304 }
305 fprintf(f, "\n");
306
307 double PctIdBestP = max(Hit.PctIdQA, Hit.PctIdQB);
308 double Div = (Hit.PctIdQM - PctIdBestP)*100.0/PctIdBestP;
309
310 unsigned LTot = Hit.CS_LY + Hit.CS_LN + Hit.CS_LA;
311 unsigned RTot = Hit.CS_RY + Hit.CS_RN + Hit.CS_RA;
312
313 double PctL = Pct(Hit.CS_LY, LTot);
314 double PctR = Pct(Hit.CS_RY, RTot);
315
316 fprintf(f,
317 "Ids. QA %.1f%%, QB %.1f%%, AB %.1f%%, QModel %.1f%%, Div. %+.1f%%\n",
318 Hit.PctIdQA,
319 Hit.PctIdQB,
320 Hit.PctIdAB,
321 Hit.PctIdQM,
322 Div);
323
324 fprintf(f,
325 "Diffs Left %u: N %u, A %u, Y %u (%.1f%%); Right %u: N %u, A %u, Y %u (%.1f%%), Score %.4f\n",
326 LTot, Hit.CS_LN, Hit.CS_LA, Hit.CS_LY, PctL,
327 RTot, Hit.CS_RN, Hit.CS_RA, Hit.CS_RY, PctR,
328 Hit.Score);
329 }
330