bm->theta = 0.85;
bm->n_hap = 2;
bm->eta = 0.03;
+ bm->cap_mapQ = 60;
return bm;
}
rms = 0;
for (j = n - 1; j >= 0; --j) { // calculate esum and fsum
uint32_t info = bm->aux->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) {
if (w[k] < 0xff) ++w[k];
++b->c[k&3];
}
- rms += (int)(info&0x7f) * (info&0x7f);
+ tmp = (int)(info&0x7f) < bm->cap_mapQ? (int)(info&0x7f) : bm->cap_mapQ;
+ rms += tmp * tmp;
}
b->rms_mapQ = (uint8_t)(sqrt((double)rms / n) + .499);
// rescale ->c[]
if (p[j<<2|k] < 0.0) p[j<<2|k] = 0.0;
}
+ { // fix p[k<<2|k]
+ float max1, max2, min1, min2;
+ int max_k, min_k;
+ max_k = min_k = -1;
+ max1 = max2 = -1.0; min1 = min2 = 1e30;
+ for (k = 0; k < 4; ++k) {
+ if (b->esum[k] > max1) {
+ max2 = max1; max1 = b->esum[k]; max_k = k;
+ } else if (b->esum[k] > max2) max2 = b->esum[k];
+ }
+ for (k = 0; k < 4; ++k) {
+ if (p[k<<2|k] < min1) {
+ min2 = min1; min1 = p[k<<2|k]; min_k = k;
+ } else if (p[k<<2|k] < min2) min2 = p[k<<2|k];
+ }
+ if (max1 > max2 && (min_k != max_k || min1 + 1.0 > min2))
+ p[max_k<<2|max_k] = min1 > 1.0? min1 - 1.0 : 0.0;
+ }
+
// convert necessary information to glf1_t
g->ref_base = ref_base; g->max_mapQ = b->rms_mapQ;
g->depth = n > 16777215? 16777215 : n;
if (types[i] != 0) { // then try the other way to calculate the score
for (ps = s = 0, l = seg.qbeg; c->pos + l + types[i] < right && l < seg.qend; ++l) {
int cq = bam1_seqi(bam1_seq(p->b), l), ct;
- ct = c->pos + l + types[i] >= left? ref2[c->pos + l + types[i] - left] : 15;
+ ct = c->pos + l - seg.qbeg + types[i] >= left? ref2[c->pos + l - seg.qbeg + types[i] - left] : 15;
if (cq < 15 && ct < 15) {
s += cq == ct? 1 : -mi->mm_penalty;
if (cq != ct) ps += bam1_qual(p->b)[l];
}
if (score[i*n+j] < s) score[i*n+j] = s; // choose the higher of the two scores
if (pscore[i*n+j] > ps) pscore[i*n+j] = ps;
- if (types[i] != 0) score[i*n+j] -= mi->indel_err;
+ //if (types[i] != 0) score[i*n+j] -= mi->indel_err;
//printf("%d, %d, %d, %d, %d, %d, %d\n", p->b->core.pos + 1, seg.qbeg, i, types[i], j,
// score[i*n+j], pscore[i*n+j]);
}
if (s1 > s2) ret->gl[0] += s1 - s2 < mi->q_indel? s1 - s2 : mi->q_indel;
else ret->gl[1] += s2 - s1 < mi->q_indel? s2 - s1 : mi->q_indel;
}
+ // write cnt_ref and cnt_ambi
+ if (max1_i != 0 && max2_i != 0) {
+ for (j = 0; j < n; ++j) {
+ int diff1 = score[j] - score[max1_i * n + j];
+ int diff2 = score[j] - score[max2_i * n + j];
+ if (diff1 > 0 && diff2 > 0) ++ret->cnt_ref;
+ else if (diff1 == 0 || diff2 == 0) ++ret->cnt_ambi;
+ }
+ }
}
free(score); free(pscore); free(ref2); free(inscns);
}