- q[0] = q[1] = r[0] = r[1] = kk[0] = kk[1] = kl[0] = kl[1] = -1;
- // find the right boundary
- for (k = 0, score = max = 0, x = c->pos, y = 0; k < c->n_cigar; ++k) {
- int op = cigar[k]&0xf;
- int ol = cigar[k]>>4;
- if (op == BAM_CMATCH) {
- for (j = 0; j < ol; ++j) {
- int c1, c2, z = y + j;
- c1 = bam_nt16_nt4_table[bam1_seqi(seq, z)];
- if (ref[x+j] == 0) return -1;
- c2 = bam_nt16_nt4_table[(int)bam_nt16_table[(int)ref[x+j]]];
- if (c1 < 3 && c2 < 3)
- score += c1 == c2? 5 : -4;
- if (score < 0) score = 0;
- if (score > max) max = score, q[1] = z, r[1] = x+j, kk[1] = k, kl[1] = j + 1;
- }
- x += ol; y += ol;
- } else if (op == BAM_CINS) {
- score -= 4 - ol * 3;
- y += ol;
- if (score < 0) score = 0;
- } else if (op == BAM_CDEL) {
- score -= 4 - ol * 3;
- x += ol;
- if (score < 0) score = 0;
- } else if (op == BAM_CSOFT_CLIP) y += ol;
- else if (op == BAM_CREF_SKIP) x += ol;
- }
- if (q[1] < 0) return -1; // no high scoring segments
- endx = x - 1;
- // find the left boundary
- for (k = c->n_cigar - 1, score = max = 0, x = x-1, y = y-1; k >= 0; --k) {
- int op = cigar[k]&0xf;
- int ol = cigar[k]>>4;
- if (op == BAM_CMATCH) {
- for (j = 0; j < ol; ++j) {
- int c1, c2, z = y - j;
- c1 = bam_nt16_nt4_table[bam1_seqi(seq, z)];
- if (ref[x+j] == 0) return -1;
- c2 = bam_nt16_nt4_table[(int)bam_nt16_table[(int)ref[x-j]]];
- if (c1 < 3 && c2 < 3)
- score += c1 == c2? 5 : -4;
- if (score < 0) score = 0;
- if (score > max) max = score, q[0] = z, r[0] = x-j, kk[0] = k, kl[0] = j + 1;
- }
- x -= ol; y -= ol;
- } else if (op == BAM_CINS) {
- score -= 4 - ol * 3;
- y -= ol;
- if (score < 0) score = 0;
- } else if (op == BAM_CDEL) {
- score -= 4 - ol * 3;
- x -= ol;
- if (score < 0) score = 0;
- } else if (op == BAM_CSOFT_CLIP) y -= ol;
- else if (op == BAM_CREF_SKIP) x -= ol;
+ kpa_par_t conf = kpa_par_def;
+ uint8_t *bq = 0, *zq = 0, *qual = bam1_qual(b);
+ if ((c->flag & BAM_FUNMAP) || b->core.l_qseq == 0) return -1; // do nothing
+ // test if BQ or ZQ is present
+ if ((bq = bam_aux_get(b, "BQ")) != 0) ++bq;
+ if ((zq = bam_aux_get(b, "ZQ")) != 0 && *zq == 'Z') ++zq;
+ if (bq && redo_baq)
+ {
+ bam_aux_del(b, bq-1);
+ bq = 0;