- { // 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);
- }