+ /* The following block fixes a long-existing flaw in the INDEL
+ * calling model: the interference of nearby SNPs. However, it also
+ * reduces the power because sometimes, substitutions caused by
+ * indels are not distinguishable from true mutations. Multiple
+ * sequence realignment helps to increase the power.
+ */
+ { // construct per-sample consensus
+ int L = right - left + 1, max_i, max2_i;
+ uint32_t *cns, max, max2;
+ char *ref0, *r;
+ ref_sample = calloc(n, sizeof(void*));
+ cns = calloc(L, 4);
+ ref0 = calloc(L, 1);
+ for (i = 0; i < right - left; ++i)
+ ref0[i] = bam_nt16_table[(int)ref[i+left]];
+ for (s = 0; s < n; ++s) {
+ r = ref_sample[s] = calloc(L, 1);
+ memset(cns, 0, sizeof(int) * L);
+ // collect ref and non-ref counts
+ for (i = 0; i < n_plp[s]; ++i) {
+ bam_pileup1_t *p = plp[s] + i;
+ bam1_t *b = p->b;
+ uint32_t *cigar = bam1_cigar(b);
+ uint8_t *seq = bam1_seq(b);
+ int x = b->core.pos, y = 0;
+ for (k = 0; k < b->core.n_cigar; ++k) {
+ int op = cigar[k]&0xf;
+ int j, l = cigar[k]>>4;
+ if (op == BAM_CMATCH) {
+ for (j = 0; j < l; ++j)
+ if (x + j >= left && x + j < right)
+ cns[x+j-left] += (bam1_seqi(seq, y+j) == ref0[x+j-left])? 1 : 0x10000;
+ x += l; y += l;
+ } else if (op == BAM_CDEL || op == BAM_CREF_SKIP) x += l;
+ else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) y += l;
+ }
+ }
+ // determine the consensus
+ for (i = 0; i < right - left; ++i) r[i] = ref0[i];
+ max = max2 = 0; max_i = max2_i = -1;
+ for (i = 0; i < right - left; ++i) {
+ if (cns[i]>>16 >= max>>16) max2 = max, max2_i = max_i, max = cns[i], max_i = i;
+ else if (cns[i]>>16 >= max2>>16) max2 = cns[i], max2_i = i;
+ }
+ if ((double)(max&0xffff) / ((max&0xffff) + (max>>16)) >= 0.7) max_i = -1;
+ if ((double)(max2&0xffff) / ((max2&0xffff) + (max2>>16)) >= 0.7) max2_i = -1;
+ if (max_i >= 0) r[max_i] = 15;
+ if (max2_i >= 0) r[max2_i] = 15;
+// for (i = 0; i < right - left; ++i) fputc("=ACMGRSVTWYHKDBN"[(int)r[i]], stderr); fputc('\n', stderr);
+ }
+ free(ref0); free(cns);
+ }