Merge branch 'develop' of github.com:samtools/bcftools into develop

Author: pd3

NEWS

@@ -23,6 +23,8 @@ Changes affecting specific commands:
 
     - Fix a bug which prevented reading fasta files containing empty lines in their entirety (#2424)
 
+    - Fix a bug which causes `--absent` miss some absent positions
+
 * bcftools csq
 
     - Add support for complex substitutions, such as AC>TAA

README.md

@@ -41,3 +41,9 @@ Please cite this paper when using BCFtools for your publications. http://samtool
     eprint = {https://academic.oup.com/gigascience/article-pdf/10/2/giab008/36332246/giab008.pdf},
 }
 ```
+
+### Support
+
+If you have found a bug or would like a new feature, please report the same in the GitHub [BCFtools](https://github.com/samtools/bcftools/issues) issue tracker.
+
+For any security related issue, please send a mail to [[email protected]](mailto:[email protected]) instead of reporting in the GitHub issue tracker.

consensus.c

@@ -430,7 +430,7 @@ static void flush_fa_buffer(args_t *args, int len)
 {
     if ( !args->fa_buf.l ) return;
     int nwr = 0;
-    while ( nwr + 60 <= args->fa_buf.l )
+    while ( nwr + 60 <= args->fa_buf.l && (!len || args->fa_ori_pos + nwr + 60 < args->fa_frz_pos) )
     {
         if ( fwrite(args->fa_buf.s+nwr,1,60,args->fp_out) != 60 ) error("Could not write: %s\n", args->output_fname);
         if ( fwrite("\n",1,1,args->fp_out) != 1 ) error("Could not write: %s\n", args->output_fname);
@@ -462,7 +462,9 @@ static void flush_fa_buffer(args_t *args, int len)
 }
 static void apply_absent(args_t *args, hts_pos_t pos)
 {
-    if ( !args->fa_buf.l || pos <= args->fa_frz_pos + 1 || pos <= args->fa_ori_pos ) return;
+    if ( !args->fa_buf.l ) return;
+    if ( pos <= args->fa_frz_pos + 1 ) return;  // if pos==frz+1, then there is no gap, ie=ib, nothing to fill
+    if ( pos <= args->fa_ori_pos ) return;
 
     int ie = pos && pos - args->fa_ori_pos + args->fa_mod_off < args->fa_buf.l ? pos - args->fa_ori_pos + args->fa_mod_off : args->fa_buf.l;
     int ib = args->fa_frz_mod < 0 ? 0 : args->fa_frz_mod;
@@ -585,7 +587,6 @@ static int iupac_set_allele(args_t *args, bcf1_t *rec)
 static void apply_variant(args_t *args, bcf1_t *rec)
 {
     static int warned_haplotype = 0;
-
     if ( args->absent_allele ) apply_absent(args, rec->pos);
     if ( rec->n_allele==1 && !args->missing_allele && !args->absent_allele ) { return; }
 

csq.c

@@ -3182,9 +3182,9 @@ int test_cds(args_t *args, bcf1_t *rec, vbuf_t *vbuf)
         for (ismpl=0; ismpl<args->smpl->n; ismpl++)
         {
             int32_t *gt = args->gt_arr + args->smpl->idx[ismpl]*ngts;
-            if ( gt[0]==bcf_gt_missing ) continue;
+            if ( bcf_gt_is_missing(gt[0]) ) continue;
 
-            if ( ngts>1 && gt[1]!=bcf_gt_missing && gt[1]!=bcf_int32_vector_end && bcf_gt_allele(gt[0])!=bcf_gt_allele(gt[1]) )
+            if ( ngts>1 && !bcf_gt_is_missing(gt[1]) && gt[1]!=bcf_int32_vector_end && bcf_gt_allele(gt[0])!=bcf_gt_allele(gt[1]) )
             {
                 if ( args->phase==PHASE_MERGE )
                 {
@@ -3206,7 +3206,7 @@ int test_cds(args_t *args, bcf1_t *rec, vbuf_t *vbuf)
 
             for (ihap=0; ihap<ngts; ihap++)
             {
-                if ( gt[ihap]==bcf_gt_missing || gt[ihap]==bcf_int32_vector_end ) continue;
+                if ( bcf_gt_is_missing(gt[ihap]) || gt[ihap]==bcf_int32_vector_end ) continue;
 
                 i = 2*ismpl + ihap;
 
@@ -3311,7 +3311,7 @@ void csq_stage(args_t *args, csq_t *csq, bcf1_t *rec)
             int32_t *gt = args->gt_arr + args->smpl->idx[i]*ngt;
             for (j=0; j<ngt; j++)
             {
-                if ( gt[j]==bcf_gt_missing || gt[j]==bcf_int32_vector_end ) continue;
+                if ( bcf_gt_is_missing(gt[j]) || gt[j]==bcf_int32_vector_end ) continue;
                 int ial = bcf_gt_allele(gt[j]);
                 if ( !ial || ial!=csq->type.vcf_ial ) continue;
                 csq_print_text(args, csq, args->smpl->idx[i],j+1);
@@ -3326,7 +3326,7 @@ void csq_stage(args_t *args, csq_t *csq, bcf1_t *rec)
         int32_t *gt = args->gt_arr + args->smpl->idx[i]*ngt;
         for (j=0; j<ngt; j++)
         {
-            if ( gt[j]==bcf_gt_missing || gt[j]==bcf_int32_vector_end ) continue;
+            if ( bcf_gt_is_missing(gt[j]) || gt[j]==bcf_int32_vector_end ) continue;
             int ial = bcf_gt_allele(gt[j]);
             if ( !ial || ial!=csq->type.vcf_ial ) continue;
 

doc/bcftools.txt

@@ -356,14 +356,17 @@ Add or remove annotations.
 
 *-c, --columns* 'list'::
     Comma-separated list of columns or tags to carry over from the annotation file
-    (see also *-a, --annotations*). If the annotation file is not a VCF/BCF,
+    (see also *-a, --annotations*).
+    {nbsp} +
+    {nbsp} +
+    *If the annotation file is not a VCF/BCF*,
     'list' describes the columns of the annotation file and must include CHROM,
     POS (or, alternatively, FROM,TO or BEG,END), and optionally REF and ALT. Unused
     columns which should be ignored can be indicated by "-".
     {nbsp} +
     {nbsp} +
-    If the annotation file is a VCF/BCF, only the edited columns/tags must be present and their
-    order does not matter. The columns ID, QUAL, FILTER, INFO and FORMAT
+    *If the annotation file is a VCF/BCF*, only the edited columns/tags must be present
+    (i.e. leave out CHROM, POS, REF, ALT) and their order does not matter. The columns ID, QUAL, FILTER, INFO and FORMAT
     can be edited, where INFO tags can be written both as "INFO/TAG" or simply "TAG",
     and FORMAT tags can be written as "FORMAT/TAG" or "FMT/TAG".
     The imported VCF annotations can be renamed as "DST_TAG:=SRC_TAG" or "FMT/DST_TAG:=FMT/SRC_TAG".

filter.c

@@ -1368,7 +1368,6 @@ static void filters_set_alt_string(filter_t *flt, bcf1_t *line, token_t *tok)
 }
 static void filters_set_nmissing(filter_t *flt, bcf1_t *line, token_t *tok)
 {
-    if ( !tok->nsamples ) error("The function %s works with FORMAT fields\n", tok->tag);
     assert(tok->usmpl);
 
     bcf_unpack(line, BCF_UN_FMT);
@@ -1400,7 +1399,7 @@ static void filters_set_nmissing(filter_t *flt, bcf1_t *line, token_t *tok)
             { \
                 type_t val = convert(&ptr[j * sizeof(type_t)]); \
                 if ( val==is_vector_end ) break; \
-                if ( val==bcf_gt_missing ) { nmissing++; break; } \
+                if ( bcf_gt_is_missing(val) ) { nmissing++; break; } \
             } \
         } \
     }
@@ -1411,7 +1410,8 @@ static void filters_set_nmissing(filter_t *flt, bcf1_t *line, token_t *tok)
         default: fprintf(stderr,"todo: type %d\n", fmt->type); exit(1); break;
     }
     #undef BRANCH
-    tok->nvalues = 1;
+    tok->nsamples = 0;
+    tok->nvalues  = 1;
     tok->values[0] = tok->tag[0]=='N' ? nmissing : (nsmpl ? (double)nmissing / nsmpl : 0);
 }
 static int func_npass(filter_t *flt, bcf1_t *line, token_t *rtok, token_t **stack, int nstack)

plugins/check-sparsity.c

@@ -206,7 +206,7 @@ static void test_region(args_t *args, char *reg)
             int8_t *ptr = (int8_t*) (fmt_gt->p + args->smpl[i]*fmt_gt->size);
             int ial = 0;
             for (ial=0; ial<fmt_gt->n; ial++)
-                if ( ptr[ial]==bcf_gt_missing || ptr[ial]==bcf_int8_vector_end ) break;
+                if ( bcf_gt_is_missing(ptr[ial]) || ptr[ial]==bcf_int8_vector_end ) break;
             if ( ial==0 ) continue;     // missing
             if ( ++args->nsites[i] < args->min_sites ) continue;
             if ( i+1<args->nsmpl )

plugins/guess-ploidy.c

@@ -167,7 +167,7 @@ void process_region_guess(args_t *args)
                 int32_t *ptr = args->arr + ismpl*ngt;
                 double *tmp = args->tmpf + ismpl*3;
 
-                if ( ptr[0]==bcf_gt_missing )
+                if ( bcf_gt_is_missing(ptr[0]) )
                 {
                     tmp[0] = -1;
                     continue;

plugins/mendelian2.c

@@ -566,12 +566,12 @@ static int parse_gt(int32_t *gt, int ngt, uint64_t *a, uint64_t *b)
 {
     *a = *b = 0;
 
-    if ( gt[0]==bcf_gt_missing || gt[0]==bcf_int32_vector_end ) return 0;
+    if ( bcf_gt_is_missing(gt[0]) || gt[0]==bcf_int32_vector_end ) return 0;
     *a |= 1<<bcf_gt_allele(gt[0]);
 
     if ( ngt==1 || gt[1]==bcf_int32_vector_end ) return 1;
 
-    if ( gt[1]==bcf_gt_missing ) return 0;
+    if ( bcf_gt_is_missing(gt[1]) ) return 0;
     *b |= 1<<bcf_gt_allele(gt[1]);
 
     return 2;

plugins/missing2ref.c

@@ -121,7 +121,7 @@ bcf1_t *process(bcf1_t *rec)
     // replace gts
     for (i=0; i<ngts; i++)
     {
-        if ( gts[i]==bcf_gt_missing )
+        if ( bcf_gt_is_missing(gts[i]) )
         {
             gts[i] = new_gt;
             changed++;