+int bcf_call_combine_gap(int n, const bcf_callret1_t *calls, bcf_call_t *call)
+{
+ int i, j, n_ins, n_del;
+ // combine annotations
+ call->ori_ref = 4;
+ memset(call->anno, 0, 16 * sizeof(int));
+ call->ins_len = call->del_len = 0; call->indelreg = 10000;
+ for (i = call->depth = 0, n_ins = n_del = 0; i < n; ++i) {
+ const bcf_callret1_t *r = calls + i;
+ if (r->depth > 0) {
+ call->depth += r->depth;
+ if (r->ins_len > 0) {
+ call->ins_len += r->ins_len;
+ ++n_ins;
+ }
+ if (r->del_len > 0) {
+ call->del_len += r->del_len;
+ ++n_del;
+ }
+ if (r->indelreg > 0 && call->indelreg > r->indelreg)
+ call->indelreg = r->indelreg;
+ for (j = 0; j < 16; ++j) call->anno[j] += r->anno[j];
+ }
+ }
+ if (call->depth == 0) return 0; // no indels
+ call->ins_len = n_ins? call->ins_len / n_ins : 0;
+ call->del_len = n_del? call->del_len / n_del : 0;
+ //
+ for (i = 0; i < 5; ++i) call->a[i] = -1;
+ call->a[0] = 0; call->a[1] = 1;
+ call->unseen = -1;
+ call->n_alleles = 2;
+ // set the PL array
+ if (call->n < n) {
+ call->n = n;
+ call->PL = realloc(call->PL, 15 * n);
+ }
+ {
+ int g[3];
+ double sum_min = 0.;
+ g[0] = 0; g[1] = 1; g[2] = 3;
+ for (i = 0; i < n; ++i) {
+ uint8_t *PL = call->PL + 3 * i;
+ const bcf_callret1_t *r = calls + i;
+ float min = 1e37;
+ for (j = 0; j < 3; ++j)
+ if (min > r->p[g[j]]) min = r->p[g[j]];
+ sum_min += min;
+ for (j = 0; j < 3; ++j) {
+ int y;
+ y = (int)(r->p[g[j]] - min + .499);
+ if (y > 255) y = 255;
+ PL[j] = y;
+ }
+ }
+ call->shift = (int)(sum_min + .499);
+ }
+ return 0;
+}
+