* put n_smpl to "bcf1_t" to simplify API a little

[samtools.git] / bcftools / call1.c

diff --git a/bcftools/call1.c b/bcftools/call1.c
index c0c24a6932c07bede41b8a3c0daba29dc82c07e7..5c6e477400845c743cb76016784297ec36e780d0 100644 (file)
--- a/bcftools/call1.c
+++ b/bcftools/call1.c
@@ -23,6 +23,7 @@ KSTREAM_INIT(gzFile, gzread, 16384)
 #define VC_HWE     128
 #define VC_KEEPALT 256
 #define VC_ACGT_ONLY 512
+#define VC_QCALL   1024
 
 typedef struct {
        int flag, prior_type, n1;
@@ -132,7 +133,7 @@ static int test16(bcf1_t *b, anno16_t *a)
        return test16_core(anno, a);
 }
 
-static void rm_info(int n_smpl, bcf1_t *b, const char *key)
+static void rm_info(bcf1_t *b, const char *key)
 {
        char *p, *q;
        if ((p = strstr(b->info, key)) == 0) return;
@@ -140,7 +141,7 @@ static void rm_info(int n_smpl, bcf1_t *b, const char *key)
        if (p > b->info && *(p-1) == ';') --p;
        memmove(p, q, b->l_str - (q - b->str));
        b->l_str -= q - p;
-       bcf_sync(n_smpl, b);
+       bcf_sync(b);
 }
 
 static int update_bcf1(int n_smpl, bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p1rst_t *pr, double pref, int flag)
@@ -152,7 +153,7 @@ static int update_bcf1(int n_smpl, bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p
 
        p_hwe = pr->g[0] >= 0.? test_hwe(pr->g) : 1.0; // only do HWE g[] is calculated
        test16(b, &a);
-       rm_info(n_smpl, b, "I16=");
+       rm_info(b, "I16=");
 
        memset(&s, 0, sizeof(kstring_t));
        kputc('\0', &s); kputs(b->ref, &s); kputc('\0', &s);
@@ -173,15 +174,16 @@ static int update_bcf1(int n_smpl, bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p
        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;
-       bcf_sync(n_smpl, b);
-       if (!is_var) bcf_shrink_alt(n_smpl, b, 1);
+       bcf_sync(b);
+       if (!is_var) bcf_shrink_alt(b, 1);
        else if (!(flag&VC_KEEPALT))
-               bcf_shrink_alt(n_smpl, b, pr->rank0 < 2? 2 : pr->rank0+1);
+               bcf_shrink_alt(b, pr->rank0 < 2? 2 : pr->rank0+1);
        return is_var;
 }
 
 int bcfview(int argc, char *argv[])
 {
+       extern int bcf_2qcall(bcf_hdr_t *h, bcf1_t *b);
        bcf_t *bp, *bout = 0;
        bcf1_t *b;
        int c;
@@ -196,7 +198,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;
-       while ((c = getopt(argc, argv, "N1:l:cHAGvLbSuP:t:p:")) >= 0) {
+       while ((c = getopt(argc, argv, "N1:l:cHAGvLbSuP:t:p:Q")) >= 0) {
                switch (c) {
                case '1': vc.n1 = atoi(optarg); break;
                case 'l': vc.fn_list = strdup(optarg); break;
@@ -212,6 +214,7 @@ int bcfview(int argc, char *argv[])
                case 'H': vc.flag |= VC_HWE; break;
                case 't': vc.theta = atof(optarg); break;
                case 'p': vc.pref = atof(optarg); break;
+               case 'Q': vc.flag |= VC_QCALL; break;
                case 'P':
                        if (strcmp(optarg, "full") == 0) vc.prior_type = MC_PTYPE_FULL;
                        else if (strcmp(optarg, "cond2") == 0) vc.prior_type = MC_PTYPE_COND2;
@@ -233,6 +236,7 @@ int bcfview(int argc, char *argv[])
                fprintf(stderr, "         -H        perform Hardy-Weinberg test (slower)\n");
                fprintf(stderr, "         -v        output potential variant sites only\n");
                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, "         -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 mutation rate [%.4lg]\n", vc.theta);
@@ -251,7 +255,7 @@ int bcfview(int argc, char *argv[])
        bp = vcf_open(argv[optind], moder);
        h = vcf_hdr_read(bp);
        bout = vcf_open("-", modew);
-       vcf_hdr_write(bout, h);
+       if (!(vc.flag & VC_QCALL)) vcf_hdr_write(bout, h);
        if (vc.flag & VC_CALL) {
                p1 = bcf_p1_init(h->n_smpl);
                if (vc.prior_file) {
@@ -273,7 +277,11 @@ int bcfview(int argc, char *argv[])
                        idx = bcf_idx_load(argv[optind]);
                        if (idx) {
                                uint64_t off;
-                               off = bcf_idx_query(idx, tid, begin, end);
+                               off = bcf_idx_query(idx, tid, begin);
+                               if (off == 0) {
+                                       fprintf(stderr, "[%s] no records in the query region.\n", __func__);
+                                       return 1; // FIXME: a lot of memory leaks...
+                               }
                                bgzf_seek(bp->fp, off, SEEK_SET);
                                bcf_idx_destroy(idx);
                        }
@@ -300,9 +308,13 @@ int bcfview(int argc, char *argv[])
                        if (!(l > begin && end > b->pos)) continue;
                }
                ++n_processed;
-               if (vc.flag & VC_CALL) {
+               if (vc.flag & VC_QCALL) { // output QCALL format; STOP here
+                       bcf_2qcall(h, b);
+                       continue;
+               }
+               if (vc.flag & VC_CALL) { // call variants
                        bcf_p1rst_t pr;
-                       bcf_gl2pl(h->n_smpl, b);
+                       bcf_gl2pl(b);
                        bcf_p1_cal(b, p1, &pr); // pr.g[3] is not calculated here
                        if (vc.flag&VC_HWE) bcf_p1_cal_g3(p1, pr.g);
                        if (n_processed % 100000 == 0) {
@@ -312,7 +324,7 @@ int bcfview(int argc, char *argv[])
                        if (pr.p_ref >= vc.pref && (vc.flag & VC_VARONLY)) continue;
                        update_bcf1(h->n_smpl, b, p1, &pr, vc.pref, vc.flag);
                }
-               if (vc.flag & VC_NO_GENO) {
+               if (vc.flag & VC_NO_GENO) { // do not output GENO fields
                        b->n_gi = 0;
                        b->fmt[0] = '\0';
                }