- // do alignment; this takes most of computing time for indel calling
- sc = ka_global_score((uint8_t*)ref2 + tbeg - left, tend - tbeg + abs(types[t]),
- (uint8_t*)query, qend - qbeg, &ap);
- score[K*n_types + t] = -sc;
+ { // do realignment; this is the bottleneck
+ const uint8_t *qual = bam1_qual(p->b), *bq;
+ uint8_t *qq;
+ qq = calloc(qend - qbeg, 1);
+ bq = (uint8_t*)bam_aux_get(p->b, "ZQ");
+ if (bq) ++bq; // skip type
+ for (l = qbeg; l < qend; ++l) {
+ qq[l - qbeg] = bq? qual[l] + (bq[l] - 64) : qual[l];
+ if (qq[l - qbeg] > 30) qq[l - qbeg] = 30;
+ if (qq[l - qbeg] < 7) qq[l - qbeg] = 7;
+ }
+ sc = kpa_glocal((uint8_t*)ref2 + tbeg - left, tend - tbeg + abs(types[t]),
+ (uint8_t*)query, qend - qbeg, qq, &apf1, 0, 0);
+ l = (int)(100. * sc / (qend - qbeg) + .499); // used for adjusting indelQ below
+ if (l > 255) l = 255;
+ score1[K*n_types + t] = score2[K*n_types + t] = sc<<8 | l;
+ if (sc > 5) {
+ sc = kpa_glocal((uint8_t*)ref2 + tbeg - left, tend - tbeg + abs(types[t]),
+ (uint8_t*)query, qend - qbeg, qq, &apf2, 0, 0);
+ l = (int)(100. * sc / (qend - qbeg) + .499);
+ if (l > 255) l = 255;
+ score2[K*n_types + t] = sc<<8 | l;
+ }
+ free(qq);
+ }