typedef struct {
int flag, prior_type;
- char *fn_list;
+ char *fn_list, *prior_file;
double theta, pref;
} viewconf_t;
p += 4;
for (i = 0; i < 4; ++i) {
d[i] = strtol(p, &p, 10);
- if (errno == EINVAL || errno == ERANGE) return -2.;
+ if (d[i] == 0 && (errno == EINVAL || errno == ERANGE)) return -2.;
++p;
}
kt_fisher_exact(d[0], d[1], d[2], d[3], &left, &right, &two);
static int update_bcf1(int n_smpl, bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p1rst_t *pr, double pref, int flag)
{
kstring_t s;
- int d[4], x, is_var = (pr->p_ref < pref);
+ int d[4], is_var = (pr->p_ref < pref);
double p_hwe, p_dp, p_ed, r = is_var? pr->p_ref : 1. - pr->p_ref;
p_hwe = test_hwe(pr->g);
}
kputc('\0', &s); kputc('\0', &s);
kputs(b->info, &s);
- if (p_dp >= 0. && d[2] + d[3] > 0) { // has at least one non-reference base
- if (b->info[0]) kputc(';', &s);
- ksprintf(&s, "AF1=%.3lf;AFE=%.3lf", 1.-pr->f_em, 1.-pr->f_exp);
- }
+ if (b->info[0]) kputc(';', &s);
+ ksprintf(&s, "AF1=%.3lf;AFE=%.3lf", 1.-pr->f_em, 1.-pr->f_exp);
if (p_hwe <= .2) ksprintf(&s, ";GC=%.2lf,%.2lf,%.2lf;HWE=%.3lf", pr->g[2], pr->g[1], pr->g[0], p_hwe);
if (p_dp >= 0. && p_dp <= .2) ksprintf(&s, ";TDP=%.3lf", p_dp);
if (p_ed >= 0. && p_ed <= .2) ksprintf(&s, ";TED=%.3lf", p_ed);
kputs(b->fmt, &s); kputc('\0', &s);
free(b->str);
b->m_str = s.m; b->l_str = s.l; b->str = s.s;
- x = (int)(r < 1e-100? 99 : -3.434 * log(r) + .499);
- b->qual = x > 99? 99 : x;
+ b->qual = r < 1e-100? 99 : -3.434 * log(r);
+ if (b->qual > 99) b->qual = 99;
bcf_sync(n_smpl, b);
return is_var;
}
if (strcmp(optarg, "full") == 0) vc.prior_type = MC_PTYPE_FULL;
else if (strcmp(optarg, "cond2") == 0) vc.prior_type = MC_PTYPE_COND2;
else if (strcmp(optarg, "flat") == 0) vc.prior_type = MC_PTYPE_FLAT;
+ else vc.prior_file = strdup(optarg);
break;
}
}
vcf_hdr_write(bout, h);
if (vc.flag & VC_CALL) {
p1 = bcf_p1_init(h->n_smpl);
- bcf_p1_init_prior(p1, vc.prior_type, vc.theta);
+ if (vc.prior_file) {
+ if (bcf_p1_read_prior(p1, vc.prior_file) < 0) {
+ fprintf(stderr, "[%s] fail to read the prior AFS.\n", __func__);
+ return 1;
+ }
+ } else bcf_p1_init_prior(p1, vc.prior_type, vc.theta);
}
if (vc.fn_list) hash = bcf_load_pos(vc.fn_list, h);
if (optind + 1 < argc && !(vc.flag&VC_VCFIN)) {
vcf_write(bout, h, b);
++n_processed;
}
+ if (vc.prior_file) free(vc.prior_file);
if (vc.flag & VC_CALL) bcf_p1_dump_afs(p1);
bcf_hdr_destroy(h);
bcf_destroy(b);