query: Allow intermixed FORMAT and other columns; stats: Check for mi…

…ssing genotypes to prevent out-of-bounds indexes; All tools: reflect upgraded bcf_idinfo_exists() changes
Alexdavv · Oct 3, 2013 · 3935d88 · 3935d88
1 parent 8f2bfbb
commit 3935d88
Show file tree

Hide file tree

Showing 4 changed files with 34 additions and 8 deletions.
diff --git a/filter.c b/filter.c
@@ -474,7 +474,7 @@ filter_t *filter_init(bcf_hdr_t *hdr, const char *str)
             if ( strcmp(".",out[j].key) )
             {
                 out[j].hdr_id = bcf_id2int(filter->hdr, BCF_DT_ID, out[j].key);
-                if ( out[j].hdr_id<0 || !bcf_idinfo_exists(filter->hdr,BCF_HL_FLT,out[j].hdr_id) )
+                if ( !bcf_idinfo_exists(filter->hdr,BCF_HL_FLT,out[j].hdr_id) )
                     error("The filter \"%s\" not present in the VCF header\n", out[j].key);
             }
             else

diff --git a/vcffilter.c b/vcffilter.c
@@ -59,7 +59,7 @@ static void init_data(args_t *args)
                 ksprintf(&flt_name,"Filter%d", ++i);
                 id = bcf_id2int(args->hdr,BCF_DT_ID,flt_name.s);
             }
-            while ( id>=0 && bcf_idinfo_exists(args->hdr,BCF_HL_FLT,id) );
+            while ( bcf_idinfo_exists(args->hdr,BCF_HL_FLT,id) );
         }
         ksprintf(&hdr_line,"##FILTER=<ID=%s,Description=\"Set if %s: %s\">", flt_name.s,args->filter_logic & FLT_INCLUDE ? "not true" : "true", args->filter_str);
         bcf_hdr_append(args->hdr, hdr_line.s);

diff --git a/vcfquery.c b/vcfquery.c
@@ -295,6 +295,28 @@ static fmt_t *register_tag(args_t *args, int type, char *key, int is_gtf)
     fmt->key   = key ? strdup(key) : NULL;
     fmt->is_gt_field = is_gtf;
     fmt->subscript = -1;
+
+    // Allow non-format tags, such as CHROM, INFO, etc., to appear amongst the format tags.
+    if ( key )
+    {
+        int id = bcf_id2int(args->header, BCF_DT_ID, key);
+        if ( fmt->type==T_FORMAT && !bcf_idinfo_exists(args->header,BCF_HL_FMT,id) )
+        {
+            if ( !strcmp("CHROM",key) ) { fmt->type = T_CHROM; }
+            else if ( !strcmp("POS",key) ) { fmt->type = T_POS; }
+            else if ( !strcmp("ID",key) ) { fmt->type = T_ID; }
+            else if ( !strcmp("REF",key) ) { fmt->type = T_REF; }
+            else if ( !strcmp("ALT",key) ) { fmt->type = T_ALT; }
+            else if ( !strcmp("QUAL",key) ) { fmt->type = T_QUAL; }
+            else if ( !strcmp("FILTER",key) ) { fmt->type = T_FILTER; }
+            else if ( id>=0 && bcf_idinfo_exists(args->header,BCF_HL_INFO,id) ) 
+            { 
+                fmt->type = T_INFO; 
+                fprintf(stderr,"Warning: Assuming INFO/%s\n", key);
+            }
+        }
+    }
+
     switch (fmt->type)
     {
         case T_CHROM: fmt->handler = &process_chrom; break;

diff --git a/vcfstats.c b/vcfstats.c
@@ -322,9 +322,9 @@ static void init_stats(args_t *args)
             stats->smpl_indels = (int *) calloc(args->files->n_smpl,sizeof(int));
             stats->smpl_dp     = (unsigned long int *) calloc(args->files->n_smpl,sizeof(unsigned long int));
             stats->smpl_ndp    = (int *) calloc(args->files->n_smpl,sizeof(int));
-            stats->smpl_sngl    = (int *) calloc(args->files->n_smpl,sizeof(int));
+            stats->smpl_sngl   = (int *) calloc(args->files->n_smpl,sizeof(int));
             #if HWE_STATS
-                stats->af_hwe = (int*) calloc(args->m_af*args->naf_hwe,sizeof(int));
+                stats->af_hwe  = (int*) calloc(args->m_af*args->naf_hwe,sizeof(int));
             #endif
         }
         idist_init(&stats->dp, args->dp_min,args->dp_max,args->dp_step);
@@ -361,6 +361,7 @@ static void destroy_stats(args_t *args)
             if (stats->qual_indels) free(stats->qual_indels);
         #endif
         #if HWE_STATS
+            //if ( args->files->n_smpl ) free(stats->af_hwe);
             free(stats->af_hwe);
         #endif
         free(stats->insertions);
@@ -643,10 +644,13 @@ static void do_sample_stats(args_t *args, stats_t *stats, bcf_sr_t *reader, int
     }
 
     #if HWE_STATS
-        float het_frac = (float)nhet_tot/(nhet_tot + nref_tot + nalt_tot);
-        int idx = het_frac*(args->naf_hwe - 1);
-        if ( line->n_allele>1 ) idx += args->naf_hwe*args->tmp_iaf[1]; 
-        stats->af_hwe[idx]++;
+        if ( nhet_tot + nref_tot + nalt_tot )
+        {
+            float het_frac = (float)nhet_tot/(nhet_tot + nref_tot + nalt_tot);
+            int idx = het_frac*(args->naf_hwe - 1);
+            if ( line->n_allele>1 ) idx += args->naf_hwe*args->tmp_iaf[1]; 
+            stats->af_hwe[idx]++;
+        }
     #endif
 
     if ( (fmt_ptr = bcf_get_fmt(reader->header,reader->buffer[0],"DP")) )