New -*,--keep-unseen-allele option. Resolves #2015

Author: pd3

NEWS

@@ -21,7 +21,14 @@ Changes affecting specific commands:
     - Fix `bcftools annotate --mark-sites`, VCF sites overlapping regions in a BED file
       were not annotated (#1989)
 
-* bcftools call and mpileup
+* bcftools call
+
+    - Output MIN_DP rather than MinDP in gVCF mode
+
+    - New `-*, --keep-unseen-allele` option to output the unobserved allele <*>,
+      intended for gVCF.
+
+* bcftools mpileup
 
     - Output MIN_DP rather than MinDP in gVCF mode
 

call.h

@@ -33,7 +33,7 @@ THE SOFTWARE.  */
 #define CALL_VARONLY        (1<<1)
 #define CALL_CONSTR_TRIO    (1<<2)
 #define CALL_CONSTR_ALLELES (1<<3)
-//
+#define CALL_KEEP_UNSEEN    (1<<4)
 #define CALL_FMT_PV4        (1<<5)
 #define CALL_FMT_GQ         (1<<6)
 #define CALL_FMT_GP         (1<<7)
@@ -125,8 +125,7 @@ call_t;
 void error(const char *format, ...);
 
 /*
- *  call() - return -1 value on critical error; -2 to skip the site; or the number of non-reference
- *            alleles on success.
+ *  call() - return -1 value on critical error; -2 to skip the site; or the number of alleles on success
  */
 int mcall(call_t *call, bcf1_t *rec);    // multiallic and rare-variant calling model
 int ccall(call_t *call, bcf1_t *rec);    // the default consensus calling model

doc/bcftools.1

@@ -2,12 +2,12 @@
 .\"     Title: bcftools
 .\"    Author: [see the "AUTHOR(S)" section]
 .\" Generator: Asciidoctor 2.0.16
-.\"      Date: 2023-08-23
+.\"      Date: 2023-10-23
 .\"    Manual: \ \&
 .\"    Source: \ \&
 .\"  Language: English
 .\"
-.TH "BCFTOOLS" "1" "2023-08-23" "\ \&" "\ \&"
+.TH "BCFTOOLS" "1" "2023-10-23" "\ \&" "\ \&"
 .ie \n(.g .ds Aq \(aq
 .el       .ds Aq '
 .ss \n[.ss] 0
@@ -51,7 +51,7 @@ standard input (stdin) and outputs to the standard output (stdout). Several
 commands can thus be  combined  with  Unix pipes.
 .SS "VERSION"
 .sp
-This manual page was last updated \fB2023\-08\-23 15:48 CEST\fP and refers to bcftools git version \fB1.18\-8\-g23dcc5a+\fP.
+This manual page was last updated \fB2023\-10\-23 09:31 CEST\fP and refers to bcftools git version \fB1.18\-19\-g6374c6a+\fP.
 .SS "BCF1"
 .sp
 The obsolete BCF1 format output by versions of samtools <= 0.1.19 is \fBnot\fP
@@ -1054,6 +1054,11 @@ output all alternate alleles present in the alignments even if they do not
 appear in any of the genotypes
 .RE
 .sp
+\fB\-\fP*\fB, \-\-keep\-unseen\-allele\fP
+.RS 4
+keep the unobserved allele <*> or <NON_REF>, useful mainly for gVCF output
+.RE
+.sp
 \fB\-f, \-\-format\-fields\fP \fIlist\fP
 .RS 4
 comma\-separated list of FORMAT fields to output for each sample. Currently
@@ -2796,7 +2801,7 @@ see \fBCommon Options\fP
 .sp
 \fB\-w, \-\-write\fP \fILIST\fP
 .RS 4
-list of input files to output given as 1\-based indices. With \fB\-p\fP and no
+comma\-separated list of input files to output given as 1\-based indices. With \fB\-p\fP and no
 \fB\-w\fP, all files are written.
 .RE
 .sp
@@ -4428,6 +4433,13 @@ continue even when some samples requested via \fB\-s/\-S\fP do not exist
 learn by example, see below
 .RE
 .sp
+\fB\-F, \-\-print\-filtered\fP \fISTR\fP
+.RS 4
+by default, samples failing \fB\-i/\-e\fP filtering expressions are suppressed from output
+when FORMAT fields are queried (for example \fI%CHROM %POS [ %GT]\fP).  With \fB\-F\fP, such
+fields will be still printed but instead of their actual value, \fISTR\fP will be used.
+.RE
+.sp
 \fB\-H, \-\-print\-header\fP
 .RS 4
 print header

doc/bcftools.html

@@ -50,7 +50,7 @@ <h2 id="_description">DESCRIPTION</h2>
 <div class="sect2">
 <h3 id="_version">VERSION</h3>
 <div class="paragraph">
-<p>This manual page was last updated <strong>2023-08-23 15:48 CEST</strong> and refers to bcftools git version <strong>1.18-8-g23dcc5a+</strong>.</p>
+<p>This manual page was last updated <strong>2023-10-23 09:31 CEST</strong> and refers to bcftools git version <strong>1.18-19-g6374c6a+</strong>.</p>
 </div>
 </div>
 <div class="sect2">
@@ -830,6 +830,10 @@ <h4 id="_inputoutput_options">Input/output options:</h4>
 <p>output all alternate alleles present in the alignments even if they do not
 appear in any of the genotypes</p>
 </dd>
+<dt class="hdlist1"><strong>-</strong>*<strong>, --keep-unseen-allele</strong></dt>
+<dd>
+<p>keep the unobserved allele &lt;*&gt; or &lt;NON_REF&gt;, useful mainly for gVCF output</p>
+</dd>
 <dt class="hdlist1"><strong>-f, --format-fields</strong> <em>list</em></dt>
 <dd>
 <p>comma-separated list of FORMAT fields to output for each sample. Currently
@@ -2487,7 +2491,7 @@ <h3 id="isec">bcftools isec [<em>OPTIONS</em>]  <em>A.vcf.gz</em> <em>B.vcf.gz</
 </dd>
 <dt class="hdlist1"><strong>-w, --write</strong> <em>LIST</em></dt>
 <dd>
-<p>list of input files to output given as 1-based indices. With <strong>-p</strong> and no
+<p>comma-separated list of input files to output given as 1-based indices. With <strong>-p</strong> and no
 <strong>-w</strong>, all files are written.</p>
 </dd>
 <dt class="hdlist1"><strong>--write-index</strong></dt>
@@ -3849,6 +3853,12 @@ <h3 id="query">bcftools query [<em>OPTIONS</em>] <em>file.vcf.gz</em> [<em>file.
 <dd>
 <p>learn by example, see below</p>
 </dd>
+<dt class="hdlist1"><strong>-F, --print-filtered</strong> <em>STR</em></dt>
+<dd>
+<p>by default, samples failing <strong>-i/-e</strong> filtering expressions are suppressed from output
+when FORMAT fields are queried (for example <em>%CHROM %POS [ %GT]</em>).  With <strong>-F</strong>, such
+fields will be still printed but instead of their actual value, <em>STR</em> will be used.</p>
+</dd>
 <dt class="hdlist1"><strong>-H, --print-header</strong></dt>
 <dd>
 <p>print header</p>
@@ -5317,7 +5327,7 @@ <h2 id="_copying">COPYING</h2>
 </div>
 <div id="footer">
 <div id="footer-text">
-Last updated 2023-08-23 15:48:11 +0200
+Last updated 2023-10-23 09:31:31 +0200
 </div>
 </div>
 </body>

doc/bcftools.txt

@@ -624,6 +624,9 @@ demand. The original calling model can be invoked with the *-c* option.
     output all alternate alleles present in the alignments even if they do not
     appear in any of the genotypes
 
+*-***, --keep-unseen-allele*::
+    keep the unobserved allele <*> or <NON_REF>, useful mainly for gVCF output
+
 *-f, --format-fields* 'list'::
     comma-separated list of FORMAT fields to output for each sample. Currently
     GQ and GP fields are supported. For convenience, the fields can be given

gvcf.c

@@ -98,7 +98,6 @@ bcf1_t *gvcf_write(gvcf_t *gvcf, htsFile *fh, bcf_hdr_t *hdr, bcf1_t *rec, int i
     // encountered, or other conditions not met (block broken by a non-ref or DP too low).
     int needs_flush = can_collapse ? 0 : 1;
 
-
     // Can the record be included in a gVCF block? That is, is this a ref-only site?
     if ( rec && can_collapse )
     {

mcall.c

@@ -1558,6 +1558,11 @@ int mcall(call_t *call, bcf1_t *rec)
     call->nals_new = 0;
     for (i=0; i<nals_ori; i++)
     {
+        if ( (call->flag&CALL_KEEP_UNSEEN) && i==unseen && call->nals_new==1 )
+        {
+            call->nals_new++;
+            call->als_new |= 1<<i;
+        }
         if ( i>0 && i==unseen ) continue;
         if ( call->flag & CALL_KEEPALT ) call->als_new |= 1<<i;
         if ( call->als_new & (1<<i) ) call->nals_new++;
@@ -1669,6 +1674,6 @@ int mcall(call_t *call, bcf1_t *rec)
 
     bcf_update_info_int32(call->hdr, rec, "I16", NULL, 0);     // remove I16 tag
 
-    return call->nals_new;
+    return is_variant ? call->nals_new : 1;
 }
 

vcfcall.c

@@ -910,6 +910,7 @@ static void usage(args_t *args)
     fprintf(stderr, "\n");
     fprintf(stderr, "Input/output options:\n");
     fprintf(stderr, "   -A, --keep-alts                 Keep all possible alternate alleles at variant sites\n");
+    fprintf(stderr, "   -*, --keep-unseen-allele        Keep the unobserved allele <*> or <NON_REF>\n");
     fprintf(stderr, "   -a, --annotate LIST             Optional tags to output (lowercase allowed); '?' to list available tags\n");
     fprintf(stderr, "   -F, --prior-freqs AN,AC         Use prior allele frequencies, determined from these pre-filled tags\n");
     fprintf(stderr, "   -G, --group-samples FILE|-      Group samples by population (file with \"sample\\tgroup\") or \"-\" for single-sample calling.\n");
@@ -987,6 +988,7 @@ int main_vcfcall(int argc, char *argv[])
         {"targets-file",required_argument,NULL,'T'},
         {"threads",required_argument,NULL,9},
         {"keep-alts",no_argument,NULL,'A'},
+        {"keep-unseen-allele",no_argument,NULL,'*'},
         {"insert-missed",no_argument,NULL,'i'},
         {"skip-Ns",no_argument,NULL,'N'},            // now the new default
         {"keep-masked-refs",no_argument,NULL,'M'},
@@ -1008,7 +1010,7 @@ int main_vcfcall(int argc, char *argv[])
     };
 
     char *tmp = NULL;
-    while ((c = getopt_long(argc, argv, "h?o:O:r:R:s:S:t:T:ANMV:vcmp:C:n:P:f:a:ig:XYF:G:", loptions, NULL)) >= 0)
+    while ((c = getopt_long(argc, argv, "h?o:O:r:R:s:S:t:T:A*NMV:vcmp:C:n:P:f:a:ig:XYF:G:", loptions, NULL)) >= 0)
     {
         switch (c)
         {
@@ -1026,6 +1028,7 @@ int main_vcfcall(int argc, char *argv[])
             case 'M': args.flag &= ~CF_ACGT_ONLY; break;     // keep sites where REF is N
             case 'N': args.flag |= CF_ACGT_ONLY; break;      // omit sites where first base in REF is N (the new default)
             case 'A': args.aux.flag |= CALL_KEEPALT; break;
+            case '*': args.aux.flag |= CALL_KEEP_UNSEEN; break;
             case 'c': args.flag |= CF_CCALL; break;          // the original EM based calling method
             case 'i': args.flag |= CF_INS_MISSED; break;
             case 'v': args.aux.flag |= CALL_VARONLY; break;