- 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]);
+ f = f0 = mc_freq0(_ref0, n_plp, plp, ma, &_ref, &_alt);
+ if (f >= 0.0) {
+ double flast = f;
+ for (j = 0; j < 10; ++j) {
+ f = mc_freq_iter(flast, ma);
+ if (fabs(f - flast) < 1e-3) break;
+ flast = f;
+ }
+ }
+ printf("%s\t%d\t.\t%c\t", h->target_name[tid], pos + 1, _ref0b);