In ksprintf(), change "%lf" and "%lg" to "%f" and "%g", respectively.

[samtools.git] / bcftools / call1.c

diff --git a/bcftools/call1.c b/bcftools/call1.c
index 15f6f8bcc7095c37ea5fb589efdce630680b07b7..012c2b32b7d9c731b1c509559b66f94320ae9cbb 100644 (file)
--- a/bcftools/call1.c
+++ b/bcftools/call1.c
@@ -25,6 +25,7 @@ KSTREAM_INIT(gzFile, gzread, 16384)
 #define VC_QCALL   1024
 #define VC_CALL_GT 2048
 #define VC_ADJLD   4096
+#define VC_NO_INDEL 8192
 
 typedef struct {
        int flag, prior_type, n1;
@@ -127,7 +128,7 @@ static int test16(bcf1_t *b, anno16_t *a)
        if ((p = strstr(b->info, "I16=")) == 0) return -1;
        p += 4;
        for (i = 0; i < 16; ++i) {
-               anno[i] = strtol(p, &p, 10);
+               errno = 0; anno[i] = strtol(p, &p, 10);
                if (anno[i] == 0 && (errno == EINVAL || errno == ERANGE)) return -2;
                ++p;
        }
@@ -149,7 +150,7 @@ static int update_bcf1(int n_smpl, bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p
 {
        kstring_t s;
        int is_var = (pr->p_ref < pref);
-       double p_hwe, r = is_var? pr->p_ref : 1. - pr->p_ref;
+       double p_hwe, r = is_var? pr->p_ref : pr->p_var, fq;
        anno16_t a;
 
        p_hwe = pr->g[0] >= 0.? test_hwe(pr->g) : 1.0; // only do HWE g[] is calculated
@@ -162,20 +163,24 @@ static int update_bcf1(int n_smpl, bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p
        kputs(b->info, &s);
        if (b->info[0]) kputc(';', &s);
 //     ksprintf(&s, "AF1=%.4lg;AFE=%.4lg;CI95=%.4lg,%.4lg", 1.-pr->f_em, 1.-pr->f_exp, pr->cil, pr->cih);
-       ksprintf(&s, "AF1=%.4lg;CI95=%.4lg,%.4lg", 1.-pr->f_em, pr->cil, pr->cih);
+       ksprintf(&s, "AF1=%.4g;CI95=%.4g,%.4g", 1.-pr->f_em, pr->cil, pr->cih);
        ksprintf(&s, ";DP4=%d,%d,%d,%d;MQ=%d", a.d[0], a.d[1], a.d[2], a.d[3], a.mq);
+       fq = pr->p_ref_folded < 0.5? -4.343 * log(pr->p_ref_folded) : 4.343 * log(pr->p_var_folded);
+       if (fq < -999) fq = -999;
+       if (fq > 999) fq = 999;
+       ksprintf(&s, ";FQ=%.3g", fq);
        if (a.is_tested) {
-               if (pr->pc[0] >= 0.) ksprintf(&s, ";PC4=%lg,%lg,%lg,%lg", pr->pc[0], pr->pc[1], pr->pc[2], pr->pc[3]);
-               ksprintf(&s, ";PV4=%.2lg,%.2lg,%.2lg,%.2lg", a.p[0], a.p[1], a.p[2], a.p[3]);
+               if (pr->pc[0] >= 0.) ksprintf(&s, ";PC4=%g,%g,%g,%g", pr->pc[0], pr->pc[1], pr->pc[2], pr->pc[3]);
+               ksprintf(&s, ";PV4=%.2g,%.2g,%.2g,%.2g", a.p[0], a.p[1], a.p[2], a.p[3]);
        }
        if (pr->g[0] >= 0. && p_hwe <= .2)
-               ksprintf(&s, ";GC=%.2lf,%.2lf,%.2lf;HWE=%.3lf", pr->g[2], pr->g[1], pr->g[0], p_hwe);
+               ksprintf(&s, ";GC=%.2f,%.2f,%.2f;HWE=%.3f", pr->g[2], pr->g[1], pr->g[0], p_hwe);
        kputc('\0', &s);
        kputs(b->fmt, &s); kputc('\0', &s);
        free(b->str);
        b->m_str = s.m; b->l_str = s.l; b->str = s.s;
-       b->qual = r < 1e-100? 99 : -4.343 * log(r);
-       if (b->qual > 99) b->qual = 99;
+       b->qual = r < 1e-100? 999 : -4.343 * log(r);
+       if (b->qual > 999) b->qual = 999;
        bcf_sync(b);
        if (!is_var) bcf_shrink_alt(b, 1);
        else if (!(flag&VC_KEEPALT))
@@ -215,7 +220,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.;
-       while ((c = getopt(argc, argv, "N1:l:cHAGvbSuP:t:p:QgLi:")) >= 0) {
+       while ((c = getopt(argc, argv, "N1:l:cHAGvbSuP:t:p:QgLi:I")) >= 0) {
                switch (c) {
                case '1': vc.n1 = atoi(optarg); break;
                case 'l': vc.fn_list = strdup(optarg); break;
@@ -229,6 +234,7 @@ int bcfview(int argc, char *argv[])
                case 'u': vc.flag |= VC_UNCOMP | VC_BCFOUT; break;
                case 'H': vc.flag |= VC_HWE; break;
                case 'g': vc.flag |= VC_CALL_GT | VC_CALL; break;
+               case 'I': vc.flag |= VC_NO_INDEL; break;
                case 't': vc.theta = atof(optarg); break;
                case 'p': vc.pref = atof(optarg); break;
                case 'i': vc.indel_frac = atof(optarg); break;
@@ -257,11 +263,12 @@ int bcfview(int argc, char *argv[])
                fprintf(stderr, "         -N        skip sites where REF is not A/C/G/T\n");
                fprintf(stderr, "         -Q        output the QCALL likelihood format\n");
                fprintf(stderr, "         -L        calculate LD for adjacent sites\n");
+               fprintf(stderr, "         -I        skip indels\n");
                fprintf(stderr, "         -1 INT    number of group-1 samples [0]\n");
                fprintf(stderr, "         -l FILE   list of sites to output [all sites]\n");
-               fprintf(stderr, "         -t FLOAT  scaled substitution mutation rate [%.4lg]\n", vc.theta);
-               fprintf(stderr, "         -i FLOAT  indel-to-substitution ratio [%.4lg]\n", vc.indel_frac);
-               fprintf(stderr, "         -p FLOAT  variant if P(ref|D)<FLOAT [%.3lg]\n", vc.pref);
+               fprintf(stderr, "         -t FLOAT  scaled substitution mutation rate [%.4g]\n", vc.theta);
+               fprintf(stderr, "         -i FLOAT  indel-to-substitution ratio [%.4g]\n", vc.indel_frac);
+               fprintf(stderr, "         -p FLOAT  variant if P(ref|D)<FLOAT [%.3g]\n", vc.pref);
                fprintf(stderr, "         -P STR    type of prior: full, cond2, flat [full]\n");
                fprintf(stderr, "\n");
                return 1;
@@ -290,7 +297,7 @@ int bcfview(int argc, char *argv[])
                        bcf_p1_set_n1(p1, vc.n1);
                        bcf_p1_init_subprior(p1, vc.prior_type, vc.theta);
                }
-               if (vc.indel_frac > 0.) bcf_p1_indel_prior(p1, vc.indel_frac);
+               if (vc.indel_frac > 0.) bcf_p1_indel_prior(p1, vc.indel_frac); // otherwise use the default indel_frac
        }
        if (vc.fn_list) hash = bcf_load_pos(vc.fn_list, h);
        if (optind + 1 < argc && !(vc.flag&VC_VCFIN)) {
@@ -311,7 +318,9 @@ int bcfview(int argc, char *argv[])
                }
        }
        while (vcf_read(bp, h, b) > 0) {
-               if (vc.flag & VC_ACGT_ONLY) {
+               int is_indel = bcf_is_indel(b);
+               if ((vc.flag & VC_NO_INDEL) && is_indel) continue;
+               if ((vc.flag & VC_ACGT_ONLY) && !is_indel) {
                        int x;
                        if (b->ref[0] == 0 || b->ref[1] != 0) continue;
                        x = toupper(b->ref[0]);
@@ -353,7 +362,7 @@ int bcfview(int argc, char *argv[])
                                kstring_t s;
                                s.m = s.l = 0; s.s = 0;
                                if (*b->info) kputc(';', &s);
-                               ksprintf(&s, "NEIR=%.3lf;NEIF=%.3lf,%.3lf", r2, f[0]+f[2], f[0]+f[1]);
+                               ksprintf(&s, "NEIR=%.3f;NEIF=%.3f,%.3f", r2, f[0]+f[2], f[0]+f[1]);
                                bcf_append_info(b, s.s, s.l);
                                free(s.s);
                        }