if (theta <= 0.) theta = CALL_DEFTHETA;
bca = calloc(1, sizeof(bcf_callaux_t));
bca->capQ = 60;
- bca->openQ = 40; bca->extQ = 20; bca->tandemQ = 80;
+ bca->openQ = 40; bca->extQ = 20; bca->tandemQ = 100;
bca->min_baseQ = min_baseQ;
bca->e = errmod_init(1. - theta);
return bca;
memset(r, 0, sizeof(bcf_callret1_t));
for (i = n = 0; i < _n; ++i) {
const bam_pileup1_t *p = pl + i;
- int q, b, mapQ, baseQ, is_diff, min_dist;
+ int q, b, mapQ, baseQ, is_diff, min_dist, seqQ;
// set base
if (p->is_del || p->is_refskip || (p->b->core.flag&BAM_FUNMAP)) continue;
baseQ = q = is_indel? p->aux&0xff : (int)bam1_qual(p->b)[p->qpos]; // base/indel quality
+ seqQ = is_indel? (p->aux>>8&0xff) : 99;
if (q < bca->min_baseQ) continue;
+ if (q > seqQ) q = seqQ;
mapQ = p->b->core.qual < bca->capQ? p->b->core.qual : bca->capQ;
if (q > mapQ) q = mapQ;
if (q > 63) q = 63;
b = bam_nt16_nt4_table[b? b : ref_base]; // b is the 2-bit base
is_diff = (ref4 < 4 && b == ref4)? 0 : 1;
} else {
- b = p->aux>>8&0x3f;
+ b = p->aux>>16&0x3f;
is_diff = (b != 0);
}
bca->bases[n++] = q<<5 | (int)bam1_strand(p->b)<<4 | b;
if (((ref4 < 4 && j < 4) || (ref4 == 4 && j < 5)) && i >= 0)
call->unseen = j, call->a[j++] = qsum[i]&3;
call->n_alleles = j;
- } else call->n_alleles = j;
+ } else {
+ call->n_alleles = j;
+ if (call->n_alleles == 1) return -1; // no reliable supporting read. stop doing anything
+ }
// set the PL array
if (call->n < n) {
call->n = n;
kputc(ref[pos], &s);
for (i = 1; i < 4; ++i) {
if (bc->a[i] < 0) break;
- if (i > 1) kputc(',', &s);
+ if (i > 1) {
+ kputc(',', &s); kputc(ref[pos], &s);
+ }
if (bca->indel_types[bc->a[i]] < 0) { // deletion
for (j = -bca->indel_types[bc->a[i]]; j < bca->indelreg; ++j)
kputc(ref[pos+1+j], &s);