- x |= (uint32_t)bam1_strand(p->b) << 18 | q << 8 | p->b->core.qual;
- if (p->b->core.qual < q) q = p->b->core.qual; // cap the overall quality at mapping quality
- x |= q << 24;
- qq = bam1_seqi(bam1_seq(p->b), p->qpos); // base
- q = bam_nt16_nt4_table[qq? qq : ref_base]; // q is the 2-bit base
- if (q < 4) x |= 1 << 21 | q << 16;
-
- bca->info[n++] = x;
- }
- ks_introsort_uint32_t(n, bca->info);
- r->depth = n;
- // generate esum and fsum
- memset(&aux, 0, sizeof(auxaux_t));
- for (j = n - 1, r->sum_Q2 = 0; j >= 0; --j) { // calculate esum and fsum
- uint32_t info = bca->info[j];
- int tmp;
- if (info>>24 < 4 && (info>>8&0x3f) != 0) info = 4<<24 | (info&0xffffff);
- k = info>>16&7;
- if (info>>24 > 0) {
- aux.esum[k&3] += bca->fk[aux.w[k]] * (info>>24);
- aux.fsum[k&3] += bca->fk[aux.w[k]];
- if (aux.w[k] + 1 < CALL_MAX) ++aux.w[k];
- ++aux.c[k&3];
- }
- tmp = (int)(info&0xff) < bca->capQ? (int)(info&0xff) : bca->capQ;
- r->sum_Q2 += tmp * tmp;
- }
- memcpy(r->esum, aux.esum, 5 * sizeof(float));
- // rescale ->c[]
- for (j = c = 0; j != 4; ++j) c += aux.c[j];
- if (c > 255) {
- for (j = 0; j != 4; ++j) aux.c[j] = (int)(254.0 * aux.c[j] / c + 0.499);
- for (j = c = 0; j != 4; ++j) c += aux.c[j];
- }
- // generate likelihood
- for (j = 0; j != 5; ++j) {
- float tmp;
- // homozygous
- for (k = 0, tmp = 0.0; k != 5; ++k)
- if (j != k) tmp += aux.esum[k];
- p[j*5+j] = tmp; // anything that is not j
- // heterozygous
- for (k = j + 1; k < 5; ++k) {
- for (i = 0, tmp = 0.0; i != 5; ++i)
- if (i != j && i != k) tmp += aux.esum[i];
- p[j*5+k] = p[k*5+j] = 3.01 * (aux.c[j] + aux.c[k]) + tmp;
+ if (q > seqQ) q = seqQ;
+ 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;
+ if (!is_indel) {
+ b = bam1_seqi(bam1_seq(p->b), p->qpos); // base
+ 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>>16&0x3f;
+ is_diff = (b != 0);