From: Heng Li <lh3@live.co.uk>
Date: Sun, 17 Apr 2011 23:32:01 +0000 (+0000)
Subject:  * fixed a bug in bedidx: input BED not sorted
X-Git-Url: https://git.donarmstrong.com/?p=samtools.git;a=commitdiff_plain;h=ce2a057a1bd9803e7e838e9c81551533172bcabb

 * 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 080b6fe..84a8e76 100644
--- 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 3380b30..8e57aa1 100644
--- 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 cae35eb..f5d2fd9 100644
--- 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 667a273..fc9cb29 100644
--- 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 9297831..722877d 100644
--- 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);
 		}
 	}