};
struct __bcf_p1aux_t {
- int n, M, n1, is_indel, is_folded;
+ int n, M, n1, is_indel;
double *q2p, *pdg; // pdg -> P(D|g)
double *phi, *phi_indel;
double *z, *zswap; // aux for afs
int PL_len;
};
-static void fold_array(int M, double *x)
-{
- int k;
- for (k = 0; k < M/2; ++k)
- x[k] = x[M-k] = (x[k] + x[M-k]) / 2.;
-}
-
void bcf_p1_indel_prior(bcf_p1aux_t *ma, double x)
{
int i;
return 0;
}
-void bcf_p1_set_folded(bcf_p1aux_t *p1a)
-{
- if (p1a->n1 < 0) {
- p1a->is_folded = 1;
- fold_array(p1a->M, p1a->phi);
- fold_array(p1a->M, p1a->phi_indel);
- }
-}
-
void bcf_p1_destroy(bcf_p1aux_t *ma)
{
if (ma) {
pc[1] = pc[1] == 1.? 99 : (int)(-4.343 * log(1. - pc[1]) + .499);
}
-static double mc_cal_afs(bcf_p1aux_t *ma)
+static double mc_cal_afs(bcf_p1aux_t *ma, double *p_ref_folded, double *p_var_folded)
{
int k;
- long double sum = 0.;
+ long double sum = 0., sum2;
double *phi = ma->is_indel? ma->phi_indel : ma->phi;
memset(ma->afs1, 0, sizeof(double) * (ma->M + 1));
mc_cal_y(ma);
+ // compute AFS
for (k = 0, sum = 0.; k <= ma->M; ++k)
sum += (long double)phi[k] * ma->z[k];
for (k = 0; k <= ma->M; ++k) {
ma->afs1[k] = phi[k] * ma->z[k] / sum;
if (isnan(ma->afs1[k]) || isinf(ma->afs1[k])) return -1.;
}
+ // compute folded variant probability
+ for (k = 0, sum = 0.; k <= ma->M; ++k)
+ sum += (long double)(phi[k] + phi[ma->M - k]) / 2. * ma->z[k];
+ for (k = 1, sum2 = 0.; k < ma->M; ++k)
+ sum2 += (long double)(phi[k] + phi[ma->M - k]) / 2. * ma->z[k];
+ *p_var_folded = sum2 / sum;
+ *p_ref_folded = (phi[k] + phi[ma->M - k]) / 2. * (ma->z[ma->M] + ma->z[0]) / sum;
+ // the expected frequency
for (k = 0, sum = 0.; k <= ma->M; ++k) {
ma->afs[k] += ma->afs1[k];
sum += k * ma->afs1[k];
if (b->n_alleles < 2) return -1; // FIXME: find a better solution
//
rst->rank0 = cal_pdg(b, ma);
- rst->f_exp = mc_cal_afs(ma);
- rst->p_ref = ma->is_folded? ma->afs1[ma->M] + ma->afs1[0] : ma->afs1[ma->M];
+ rst->f_exp = mc_cal_afs(ma, &rst->p_ref_folded, &rst->p_var_folded);
+ rst->p_ref = ma->afs1[ma->M];
+ for (k = 0, sum = 0.; k < ma->M; ++k)
+ sum += ma->afs1[k];
+ rst->p_var = (double)sum;
// calculate f_flat and f_em
for (k = 0, sum = 0.; k <= ma->M; ++k)
sum += (long double)ma->z[k];
void bcf_p1_dump_afs(bcf_p1aux_t *ma)
{
int k;
- if (ma->is_folded) fold_array(ma->M, ma->afs);
fprintf(stderr, "[afs]");
for (k = 0; k <= ma->M; ++k)
fprintf(stderr, " %d:%.3lf", k, ma->afs[ma->M - k]);