}
for (i = n = 0; i < _n; ++i) {
const bam_pileup1_t *p = pl + i;
- uint32_t q, x = 0;
+ uint32_t q, x = 0, qq;
if (p->is_del || (p->b->core.flag&BAM_FUNMAP)) continue;
q = (uint32_t)bam1_qual(p->b)[p->qpos];
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;
x |= q << 24;
- q = bam_nt16_nt4_table[bam1_seqi(bam1_seq(p->b), p->qpos)];
+ qq = bam1_seqi(bam1_seq(p->b), p->qpos);
+ q = bam_nt16_nt4_table[qq? qq : ref_base];
if (!p->is_del && q < 4) x |= 1 << 21 | q << 16;
bm->aux->info[n++] = x;
}
bam_segreg(pos, c, cigar, &seg);
for (ps = s = 0, l = seg.qbeg; c->pos + l < right && l < seg.qend; ++l) {
int cq = bam1_seqi(bam1_seq(p->b), l), ct;
- ct = c->pos + l >= left? ref2[c->pos + l - left] : 15; // "<" will happen if reads are too long
+ // in the following line, "<" will happen if reads are too long
+ ct = c->pos + l - seg.qbeg >= left? ref2[c->pos + l - seg.qbeg - left] : 15;
if (cq < 15 && ct < 15) {
s += cq == ct? 1 : -mi->mm_penalty;
if (cq != ct) ps += bam1_qual(p->b)[l];
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;
- //printf("%d, %d, %d, %d\n", i, types[i], j, score[i*n+j]);
+ //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]);
}
}
{ // get final result
ret->s[1][k+1] = ref[pos + k + 1];
} else ret->s[1][0] = '*';
// write count
- for (j = 0; j < n; ++j) {
- if (score[max1_i*n+j] < 0 && score[max2_i*n+j] < 0) ++ret->cnt_anti;
- else {
- int diff = score[max1_i*n+j] - score[max2_i*n+j];
- if (diff > mi->ambi_thres) ++ret->cnt1;
- else if (diff < -mi->ambi_thres) ++ret->cnt2;
- else ++ret->cnt_ambi;
- }
+ for (i = 0; i < n; ++i) {
+ const bam_pileup1_t *p = pl + i;
+ if (p->indel == ret->indel1) ++ret->cnt1;
+ else if (p->indel == ret->indel2) ++ret->cnt2;
+ else ++ret->cnt_anti;
}
// write gl[]
ret->gl[0] = ret->gl[1] = 0;
for (j = 0; j < n; ++j) {
int s1 = pscore[max1_i*n + j], s2 = pscore[max2_i*n + j];
- //printf("%d, %d\n", s1, s2);
+ //printf("%d, %d, %d, %d, %d\n", pl[j].b->core.pos+1, max1_i, max2_i, s1, s2);
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;
}