ks_tokaux_t aux;
// set ref, alt, flt, info, fmt
b->ref = b->alt = b->flt = b->info = b->fmt = 0;
- for (p = b->str, n = 0; p < b->str + b->l_str; ++p)
- if (*p == 0 && p+1 != b->str + b->l_str) tmp[n++] = p + 1;
+ for (p = b->str, n = 0; p < b->str + b->l_str; ++p) {
+ if (*p == 0 && p+1 != b->str + b->l_str) {
+ if (n == 5) {
+ ++n;
+ break;
+ } else tmp[n++] = p + 1;
+ }
+ }
if (n != 5) {
fprintf(stderr, "[%s] incorrect number of fields (%d != 5) at %d:%d\n", __func__, n, b->tid, b->pos);
return -1;
tid = begin = end = -1;
memset(&vc, 0, sizeof(viewconf_t));
- vc.prior_type = vc.n1 = -1; vc.theta = 1e-3; vc.pref = 0.5; vc.indel_frac = -1.; vc.n_perm = 0; vc.min_perm_p = 0.01; vc.min_smpl_frac = 0; vc.min_lrt = 0.01;
+ vc.prior_type = vc.n1 = -1; vc.theta = 1e-3; vc.pref = 0.5; vc.indel_frac = -1.; vc.n_perm = 0; vc.min_perm_p = 0.01; vc.min_smpl_frac = 0; vc.min_lrt = 1;
while ((c = getopt(argc, argv, "FN1:l:cC:eHAGvbSuP:t:p:QgLi:IMs:D:U:X:d:")) >= 0) {
switch (c) {
case '1': vc.n1 = atoi(optarg); break;
if (vc.flag & (VC_CALL|VC_ADJLD)) bcf_gl2pl(b);
if (vc.flag & VC_CALL) { // call variants
bcf_p1rst_t pr;
- bcf_p1_cal(b, (em[7] >= 0 && em[7] < vc.min_lrt), p1, &pr);
+ int calret = bcf_p1_cal(b, (em[7] >= 0 && em[7] < vc.min_lrt), p1, &pr);
if (n_processed % 100000 == 0) {
fprintf(stderr, "[%s] %ld sites processed.\n", __func__, (long)n_processed);
bcf_p1_dump_afs(p1);
}
pr.perm_rank = n;
}
- update_bcf1(b, p1, &pr, vc.pref, vc.flag, em);
+ if (calret >= 0) update_bcf1(b, p1, &pr, vc.pref, vc.flag, em);
} else if (vc.flag & VC_EM) update_bcf1(b, 0, 0, 0, vc.flag, em);
if (vc.flag & VC_ADJLD) { // compute LD
double f[4], r2;
minaux1_t *a = (minaux1_t*)data;
double p = 1., l = 0., f3[3];
int i;
-// printf("%lg\n", f);
+// printf("brent %lg\n", f);
if (f < 0 || f > 1) return 1e300;
f3[0] = (1.-f)*(1.-f); f3[1] = 2.*f*(1.-f); f3[2] = f*f;
for (i = a->beg; i < a->end; ++i) {
{
double f0 = *f, f3[3], err;
int i;
+// printf("em %lg\n", *f);
f3[0] = (1.-f0)*(1.-f0); f3[1] = 2.*f0*(1.-f0); f3[2] = f0*f0;
for (i = beg, f0 = 0.; i < end; ++i) {
const double *pdg = _pdg + i * 3;
if (flag & 0xf<<1) flag |= 0xf<<1;
n = b->n_smpl; n2 = n - n1;
pdg = get_pdg3(b);
+ if (pdg == 0) return -1;
for (i = 0; i < 9; ++i) x[i] = -1.;
{
if ((x[0] = est_freq(n, pdg)) < 0.) {
break;
}
}
+ if (i == b->n_gi) return -1; // no PL
if (b->n_alleles < 2) return -1; // FIXME: find a better solution
//
rst->rank0 = cal_pdg(b, ma);
#include <stdio.h>
#include <zlib.h>
+#include "ksort.h"
+KSORT_INIT_GENERIC(uint64_t)
+
#include "kseq.h"
KSTREAM_INIT(gzFile, gzread, 8192)
if (kh_exist(h, k)) {
bed_reglist_t *p = &kh_val(h, k);
if (p->idx) free(p->idx);
+ ks_introsort(uint64_t, p->n, p->a);
p->idx = bed_index_core(p->n, p->a, &p->m);
}
}