+ const char *p;
+ kputc('\t', &s);
+ if ((p = strstr(fn[i], ".bam")) != 0)
+ kputsn(fn[i], p - fn[i], &s);
+ else kputs(fn[i], &s);
+ }
+ puts(s.s);
+ free(s.s);
+ }
+ // mpileup
+ if (conf->flag & MPLP_VCF) {
+ ma = mc_init(n);
+ mc_init_prior(ma, conf->prior_type, conf->theta);
+ }
+ ref_tid = -1; ref = 0;
+ iter = bam_mplp_init(n, mplp_func, (void**)data);
+ while (bam_mplp_auto(iter, &tid, &pos, n_plp, plp) > 0) {
+ if (conf->reg && (pos < beg0 || pos >= end0)) continue; // out of the region requested
+ if (hash) {
+ khint_t k;
+ k = kh_get(64, hash, (uint64_t)tid<<32 | pos);
+ if (k == kh_end(hash)) continue;
+ }
+ if (tid != ref_tid) {
+ free(ref);
+ if (conf->fai) ref = fai_fetch(conf->fai, h->target_name[tid], &ref_len);
+ ref_tid = tid;
+ }
+ 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) {
+ putchar("ACGTN"[r.alt]);
+ if (r.alt2 >= 0 && r.alt2 < 4) printf(",%c", "ACGT"[r.alt2]);
+ } 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];