#define CALL_ETA 0.03f
#define CALL_MAX 256
#define CALL_DEFTHETA 0.83f
+#define DEF_MAPQ 20
#define CAP_DIST 25
bca->openQ = 40; bca->extQ = 20; bca->tandemQ = 100;
bca->min_baseQ = min_baseQ;
bca->e = errmod_init(1. - theta);
+ bca->min_frac = 0.002;
+ bca->min_support = 1;
return bca;
}
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;
+ mapQ = p->b->core.qual < 255? p->b->core.qual : DEF_MAPQ; // special case for mapQ==255
+ mapQ = mapQ < bca->capQ? mapQ : bca->capQ;
if (q > mapQ) q = mapQ;
if (q > 63) q = 63;
if (q < 4) q = 4;
}
bca->bases[n++] = q<<5 | (int)bam1_strand(p->b)<<4 | b;
// collect annotations
- r->qsum[b] += q;
+ if (b < 4) r->qsum[b] += q;
++r->anno[0<<2|is_diff<<1|bam1_strand(p->b)];
min_dist = p->b->core.l_qseq - 1 - p->qpos;
if (min_dist > p->qpos) min_dist = p->qpos;