- uint32_t *cigar;
- bam1_core_t *c = &p->b->core;
- int s, ps;
- bam_segreg_t seg;
- if (c->flag&BAM_FUNMAP) continue;
- cigar = bam1_cigar(p->b);
- 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;
- // 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];
- }
- }
- score[i*n + j] = s; pscore[i*n + j] = ps;
- 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 - 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];
+ int qbeg, qend, tbeg, tend;
+ if (p->b->core.flag & BAM_FUNMAP) continue;
+ qbeg = bam_tpos2qpos(&p->b->core, bam1_cigar(p->b), left, &tbeg);
+ qend = bam_tpos2qpos(&p->b->core, bam1_cigar(p->b), right, &tend);
+ for (l = qbeg; l < qend; ++l)
+ rs[l - qbeg] = bam_nt16_nt4_table[bam1_seqi(bam1_seq(p->b), l)];
+ {
+ int x, y, n_acigar, ps;
+ uint32_t *acigar;
+ ps = 0;
+ acigar = ka_global_core((uint8_t*)ref2 + tbeg - left, tend - tbeg + types[i], (uint8_t*)rs, qend - qbeg, &ap, &score[i*n+j], &n_acigar);
+ x = tbeg - left; y = 0;
+ for (l = 0; l < n_acigar; ++l) {
+ int op = acigar[l]&0xf;
+ int len = acigar[l]>>4;
+ if (op == BAM_CMATCH) {
+ int k;
+ for (k = 0; k < len; ++k)
+ if (ref2[x+k] != rs[y+k]) ps += bam1_qual(p->b)[y+k];
+ x += len; y += len;
+ } else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) {
+ if (op == BAM_CINS) ps += mi->q_indel * len;
+ y += len;
+ } else if (op == BAM_CDEL) {
+ ps += mi->q_indel * len;
+ x += len;