+ 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->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->vcf) {
+ double f0, f, fpost, pref = -1.0; // the reference allele frequency
+ int j, _ref, _alt, _ref0, depth, rms_q, _ref0b, is_var = 0, qref = 0;
+ uint64_t sqr_sum;
+ _ref0 = _ref0b = (ref && pos < ref_len)? ref[pos] : 'N';
+ _ref0 = bam_nt16_nt4_table[bam_nt16_table[_ref0]];
+ f = f0 = fpost = mc_freq0(_ref0, n_plp, plp, ma, &_ref, &_alt);
+ if (f >= 0.0) {
+ double q, flast = f;
+ for (j = 0; j < 10; ++j) {
+ f = mc_freq_iter(flast, ma);
+// printf("%lg->%lg\n", flast, f);
+ if (fabs(f - flast) < 1e-3) break;
+ flast = f;
+ }
+ pref = mc_ref_prob(ma);
+ fpost = mc_freq_post(ma);
+ if (1. - f < 1e-4) f = 1.;
+ if (1. - fpost < 1e-4) fpost = 1.;
+ is_var = (pref < .5);
+ q = is_var? pref : 1. - pref;
+ if (q < 1e-308) q = 1e-308;
+ qref = (int)(-3.434 * log(q) + .499);
+ if (qref > 99) qref = 99;
+ }
+ if (conf->var_only && !is_var) continue;
+ printf("%s\t%d\t.\t%c\t", h->target_name[tid], pos + 1, _ref0b);
+ if (is_var) {
+ if (_ref0 == _ref) putchar("ACGTN"[_alt]);
+ else printf("%c,%c", "ACGTN"[_ref], "ACGTN"[_alt]);
+ } else putchar('.');
+ printf("\t%d\t", qref);
+ if (f0 < 0.) printf("Q13\t");
+ else printf(".\t");
+ for (i = depth = 0, sqr_sum = 0; i < n; ++i) {
+ depth += n_plp[i];