* fixed a bug in bedidx: input BED not sorted

 * make bcftools a little bit more robust given wrong inputs
 * perform posterior-chi^2 by default

diff --git a/bcftools/bcf.c b/bcftools/bcf.c
index 080b6fefcfc6ba8af057bd0a9792f620ac061f7a..84a8e767ec58adb397630f609a52670779f99c89 100644 (file)
--- a/bcftools/bcf.c
+++ b/bcftools/bcf.c
@@ -103,8 +103,14 @@ int bcf_sync(bcf1_t *b)
        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;

diff --git a/bcftools/call1.c b/bcftools/call1.c
index 3380b3098ace6ff0c28fd9f7fb5951e57a7dd5b6..8e57aa1d6da3c57bd9c5dbae39385816af163773 100644 (file)
--- a/bcftools/call1.c
+++ b/bcftools/call1.c
@@ -284,7 +284,7 @@ int bcfview(int argc, char *argv[])
 
        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;
@@ -462,7 +462,7 @@ int bcfview(int argc, char *argv[])
                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);
@@ -485,7 +485,7 @@ int bcfview(int argc, char *argv[])
                                }
                                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;

diff --git a/bcftools/em.c b/bcftools/em.c
index cae35ebbec8014e11fb5b4cb5dfd9f70f77018a0..f5d2fd99e0e4d94ceb64a3d09cbc931515a49fec 100644 (file)
--- a/bcftools/em.c
+++ b/bcftools/em.c
@@ -76,7 +76,7 @@ static double prob1(double f, void *data)
        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) {
@@ -92,6 +92,7 @@ static double freq_iter(double *f, const double *_pdg, int beg, int end)
 {
        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;
@@ -178,6 +179,7 @@ int bcf_em1(const bcf1_t *b, int n1, int flag, double x[9])
        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.) {

diff --git a/bcftools/prob1.c b/bcftools/prob1.c
index 667a2732555dae4df733f12859c646eea6d63488..fc9cb29911c1abc746e63b9001a23e853aaa58f5 100644 (file)
--- a/bcftools/prob1.c
+++ b/bcftools/prob1.c
@@ -488,6 +488,7 @@ int bcf_p1_cal(const bcf1_t *b, int do_contrast, bcf_p1aux_t *ma, bcf_p1rst_t *r
                        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);

diff --git a/bedidx.c b/bedidx.c
index 92978318d4ccfbbabeeddd150458245f669062d8..722877dcc2aa1d414460c328e6e022a22befecfe 100644 (file)
--- a/bedidx.c
+++ b/bedidx.c
@@ -4,6 +4,9 @@
 #include <stdio.h>
 #include <zlib.h>
 
+#include "ksort.h"
+KSORT_INIT_GENERIC(uint64_t)
+
 #include "kseq.h"
 KSTREAM_INIT(gzFile, gzread, 8192)
 
@@ -53,6 +56,7 @@ void bed_index(void *_h)
                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);
                }
        }