- if (pos%10 == 0) mvprintw(0, tv->ccol, "%-d", pos+1);
- // print consensus
- call = bam_maqcns_call(n, pl, tv->bmc);
+ if (pos%10 == 0 && tv->mcol - tv->ccol >= 10) mvprintw(0, tv->ccol, "%-d", pos+1);
+ { // call consensus
+ bcf_callret1_t bcr;
+ int qsum[4], a1, a2, tmp;
+ double p[3], prior = 30;
+ bcf_call_glfgen(n, pl, bam_nt16_table[rb], tv->bca, &bcr);
+ for (i = 0; i < 4; ++i) qsum[i] = bcr.qsum[i]<<2 | i;
+ for (i = 1; i < 4; ++i) // insertion sort
+ for (j = i; j > 0 && qsum[j] > qsum[j-1]; --j)
+ tmp = qsum[j], qsum[j] = qsum[j-1], qsum[j-1] = tmp;
+ a1 = qsum[0]&3; a2 = qsum[1]&3;
+ p[0] = bcr.p[a1*5+a1]; p[1] = bcr.p[a1*5+a2] + prior; p[2] = bcr.p[a2*5+a2];
+ if ("ACGT"[a1] != toupper(rb)) p[0] += prior + 3;
+ if ("ACGT"[a2] != toupper(rb)) p[2] += prior + 3;
+ if (p[0] < p[1] && p[0] < p[2]) call = (1<<a1)<<16 | (int)((p[1]<p[2]?p[1]:p[2]) - p[0] + .499);
+ else if (p[2] < p[1] && p[2] < p[0]) call = (1<<a2)<<16 | (int)((p[0]<p[1]?p[0]:p[1]) - p[2] + .499);
+ else call = (1<<a1|1<<a2)<<16 | (int)((p[0]<p[2]?p[0]:p[2]) - p[1] + .499);
+ }