- // 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, %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;
+ { // write gl[]
+ int tmp, seq_err = 0;
+ double x = 1.0;
+ tmp = max1_i - max2_i;
+ if (tmp < 0) tmp = -tmp;
+ for (j = 0; j < tmp + 1; ++j) x *= INDEL_EXT_DEP;
+ seq_err = mi->q_indel * (1.0 - x) / (1.0 - INDEL_EXT_DEP);
+ 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, %d, %d, %d\n", pl[j].b->core.pos+1, max1_i, max2_i, s1, s2);
+ if (s1 > s2) ret->gl[0] += s1 - s2 < seq_err? s1 - s2 : seq_err;
+ else ret->gl[1] += s2 - s1 < seq_err? s2 - s1 : seq_err;
+ }
+ }
+ // 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;
+ }