- _ref0 = bam_nt16_nt4_table[bam_nt16_table[_ref0]];
- f = mc_freq0(_ref0, n_plp, plp, ma, &_ref, &_alt);
- printf("%s\t%d\t.\t%c\t", h->target_name[tid], pos + 1, "ACGTN"[_ref0]);
- if (_ref0 == _ref) putchar("ACGTN"[_alt]);
- else printf("%c,%c", "ACGTN"[_ref], "ACGTN"[_alt]);
- printf("\t0\tPASS\t"); // FIXME: currently these not available
- 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;
+ ref16 = bam_nt16_table[_ref0];
+ for (i = 0; i < gplp.n; ++i)
+ bcf_call_glfgen(gplp.n_plp[i], gplp.plp[i], ref16, bca, bcr + i);
+ bcf_call_combine(gplp.n, bcr, ref16, &bc);
+ bcf_call2bcf(tid, pos, &bc, b, (conf->flag&(MPLP_FMT_DP|MPLP_FMT_SP))? bcr : 0,
+ (conf->flag&MPLP_FMT_SP), 0, 0);
+ bcf_write(bp, bh, b);
+ bcf_destroy(b);
+ if (!(conf->flag&MPLP_NO_INDEL)) {
+ // call MNPs
+ if (bcf_call_mnp_prep(gplp.n, gplp.n_plp, gplp.plp, pos, bca, ref) >= 0) {
+ for (i = 0; i < gplp.n; ++i)
+ bcf_call_glfgen(gplp.n_plp[i], gplp.plp[i], B2B_REF_MNP, bca, bcr + i);
+ if (bcf_call_combine(gplp.n, bcr, B2B_REF_MNP, &bc) >= 0) {
+ b = calloc(1, sizeof(bcf1_t));
+ bcf_call2bcf(tid, pos, &bc, b, (conf->flag&(MPLP_FMT_DP|MPLP_FMT_SP))? bcr : 0,
+ (conf->flag&MPLP_FMT_SP), bca, ref);
+ bcf_write(bp, bh, b);
+ bcf_destroy(b);
+ bca->last_mnp_pos = pos;
+ }
+ }
+ // call indels
+ if (bcf_call_gap_prep(gplp.n, gplp.n_plp, gplp.plp, pos, bca, ref, rghash) >= 0) {
+ for (i = 0; i < gplp.n; ++i)
+ bcf_call_glfgen(gplp.n_plp[i], gplp.plp[i], B2B_REF_INDEL, bca, bcr + i);
+ if (bcf_call_combine(gplp.n, bcr, B2B_REF_INDEL, &bc) >= 0) {
+ b = calloc(1, sizeof(bcf1_t));
+ bcf_call2bcf(tid, pos, &bc, b, (conf->flag&(MPLP_FMT_DP|MPLP_FMT_SP))? bcr : 0,
+ (conf->flag&MPLP_FMT_SP), bca, ref);
+ bcf_write(bp, bh, b);
+ bcf_destroy(b);
+ }