- if (conf->flag & MPLP_VCF) {
- mc_rst_t r;
- int j, _ref0, depth, rms_q, _ref0b, is_var = 0, qref = 0, level = 2, tot;
- uint64_t sqr_sum;
- if (conf->flag & MPLP_AFS) level = 3;
- _ref0 = _ref0b = (ref && pos < ref_len)? ref[pos] : 'N';
- _ref0 = bam_nt16_nt4_table[bam_nt16_table[_ref0]];
- tot = mc_cal(_ref0, n_plp, plp, ma, &r, level);
- if (tot) { // has good bases
- double q;
- is_var = (r.p_ref < .5);
- q = is_var? r.p_ref : 1. - r.p_ref;
- if (q < 1e-308) q = 1e-308;
- qref = (int)(-3.434 * log(q) + .499);
- if (qref > 99) qref = 99;
- }
- if ((conf->flag & MPLP_VAR) && !is_var) continue;
- ++N; // number of processed lines
- printf("%s\t%d\t.\t%c\t", h->target_name[tid], pos + 1, _ref0b);
- if (is_var) {
- if (_ref0 == r.ref) putchar("ACGTN"[r.alt]);
- else printf("%c,%c", "ACGTN"[r.ref], "ACGTN"[r.alt]);
- } else putchar('.');
- printf("\t%d\t", qref);
- if (!tot) printf("Q13\t");
- else printf(".\t");
- for (i = depth = 0, sqr_sum = 0; i < n; ++i) {
- depth += n_plp[i];
- for (j = 0; j < n_plp[i]; ++j) {
- int q = plp[i][j].b->core.qual;
- if (q > conf->max_mq) q = conf->max_mq;
- sqr_sum += q * q;
- }
- }
- rms_q = (int)(sqrt((double)sqr_sum / depth) + .499);
- printf("DP=%d;MQ=%d", depth, rms_q);
- if (tot) {
- printf(";AF=%.3lf", 1. - r.f_em);
- if (level >= 2) printf(";AFE=%.3lf", 1-r.f_exp);
- if (conf->flag & MPLP_AFALL) {
- printf(";AF0=%.3lf;AFN=%.3lf", 1-r.f_naive, 1-r.f_nielsen);
- if (conf->flag & MPLP_AFS) printf(";AFB=%.3lf", 1-r.f_map);
- }
- }
- printf("\tGT:GQ:DP");
- if (tot) {
- for (i = 0; i < n; ++i) {
- int x = mc_call_gt(ma, r.f_exp, i);
- printf("\t%c/%c:%d:%d", "10"[((x&3)==2)], "10"[((x&3)>0)], x>>2, n_plp[i]);
- }
- } else for (i = 0; i < n; ++i) printf("\t./.:0:0");
- putchar('\n');
- if (N % MPLP_AFS_BLOCK == 0) mc_dump_afs(ma);
+ if (conf->flag & MPLP_GLF) {
+ int _ref0, ref16;
+ bcf1_t *b = calloc(1, sizeof(bcf1_t));
+ _ref0 = (ref && pos < ref_len)? ref[pos] : 'N';
+ ref16 = bam_nt16_table[_ref0];
+ for (i = 0; i < n; ++i)
+ bcf_call_glfgen(n_plp[i], plp[i], ref16, bca, bcr + i);
+ bcf_call_combine(n, bcr, ref16, &bc);
+ bcf_call2bcf(tid, pos, &bc, b);
+ bcf_write(bp, bh, b);
+ //fprintf(stderr, "%d,%d,%d\n", b->tid, b->pos, b->l_str);
+ bcf_destroy(b);