+ /* 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;
+ uint32_t *cns;
+ 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] = (cns[i] && (double)(cns[i]&0xffff) / ((cns[i]&0xffff)+(cns[i]>>16&0xffff)) < 0.7)? 15 : ref0[i];
+// for (i = 0; i < right - left; ++i) fputc("=ACMGRSVTWYHKDBN"[(int)r[i]], stderr); fputc('\n', stderr);
+ }
+ free(ref0); free(cns);
+ }