Add support for optional removal of the unseen allele

The `bcftools merge -m` command can now suppress the output of the unseen allele
<*> or <NON_REF> at variant sites (e.g. `-m both,*`) or all sites (e.g. `-m both,**`).

The `bcftools view` has a new option `-A, --trim-unseen-allele` to remove the unseen allele
at variant sites (`-A`) or all sites (`-AA`)

Resolves #2023

Author: pd3

Makefile

@@ -248,7 +248,7 @@ vcffilter.o: vcffilter.c $(htslib_vcf_h) $(htslib_synced_bcf_reader_h) $(htslib_
 vcfgtcheck.o: vcfgtcheck.c $(htslib_vcf_h) $(htslib_synced_bcf_reader_h) $(htslib_vcfutils_h) $(htslib_kbitset_h) $(htslib_hts_os_h) $(htslib_bgzf_h) $(bcftools_h) extsort.h filter.h
 vcfindex.o: vcfindex.c $(htslib_vcf_h) $(htslib_tbx_h) $(htslib_kstring_h) $(htslib_bgzf_h) $(bcftools_h)
 vcfisec.o: vcfisec.c $(htslib_vcf_h) $(htslib_synced_bcf_reader_h) $(htslib_vcfutils_h) $(htslib_hts_os_h) $(bcftools_h) $(filter_h)
-vcfmerge.o: vcfmerge.c $(htslib_vcf_h) $(htslib_synced_bcf_reader_h) $(htslib_vcfutils_h) $(htslib_faidx_h) regidx.h $(bcftools_h) vcmp.h $(htslib_khash_h)
+vcfmerge.o: vcfmerge.c $(htslib_vcf_h) $(htslib_synced_bcf_reader_h) $(htslib_vcfutils_h) $(htslib_faidx_h) regidx.h $(bcftools_h) vcmp.h $(htslib_khash_h) $(htslib_kbitset_h)
 vcfnorm.o: vcfnorm.c $(htslib_vcf_h) $(htslib_synced_bcf_reader_h) $(htslib_faidx_h) $(htslib_khash_str2int_h) $(bcftools_h) rbuf.h abuf.h gff.h regidx.h
 vcfquery.o: vcfquery.c $(htslib_vcf_h) $(htslib_synced_bcf_reader_h) $(htslib_khash_str2int_h) $(htslib_vcfutils_h) $(bcftools_h) $(filter_h) $(convert_h) $(smpl_ilist_h)
 vcfroh.o: vcfroh.c $(htslib_vcf_h) $(htslib_synced_bcf_reader_h) $(htslib_kstring_h) $(htslib_kseq_h) $(htslib_bgzf_h) $(bcftools_h) HMM.h $(smpl_ilist_h) $(filter_h)
@@ -257,7 +257,7 @@ vcfhead.o: vcfhead.c $(htslib_kstring_h) $(htslib_vcf_h) $(bcftools_h)
 vcfsom.o: vcfsom.c $(htslib_vcf_h) $(htslib_synced_bcf_reader_h) $(htslib_vcfutils_h) $(htslib_hts_os_h) $(bcftools_h)
 vcfsort.o: vcfsort.c $(htslib_vcf_h) $(htslib_kstring_h) $(htslib_hts_os_h) kheap.h $(bcftools_h)
 vcfstats.o: vcfstats.c $(htslib_vcf_h) $(htslib_synced_bcf_reader_h) $(htslib_vcfutils_h) $(htslib_faidx_h) $(bcftools_h) $(filter_h) bin.h dist.h
-vcfview.o: vcfview.c $(htslib_vcf_h) $(htslib_synced_bcf_reader_h) $(htslib_vcfutils_h) $(bcftools_h) $(filter_h) $(htslib_khash_str2int_h)
+vcfview.o: vcfview.c $(htslib_vcf_h) $(htslib_synced_bcf_reader_h) $(htslib_vcfutils_h) $(bcftools_h) $(filter_h) $(htslib_khash_str2int_h) $(htslib_kbitset_h)
 reheader.o: reheader.c $(htslib_vcf_h) $(htslib_bgzf_h) $(htslib_tbx_h) $(htslib_kseq_h) $(htslib_thread_pool_h) $(htslib_faidx_h) $(htslib_khash_str2int_h) $(bcftools_h) $(khash_str2str_h)
 tabix.o: tabix.c $(htslib_bgzf_h) $(htslib_tbx_h)
 ccall.o: ccall.c $(htslib_kfunc_h) $(call_h) kmin.h $(prob1_h)

NEWS

@@ -55,6 +55,10 @@ Changes affecting specific commands:
 
         - better description of fields
 
+* bcftools merge
+
+    - Add `-m` modifiers to suppress the output of the unseen allele <*> or <NON_REF>
+      at variant sites (e.g. `-m both,*`) or all sites (e.g. `-m both,**`)
 
 * bcftools mpileup
 
@@ -93,6 +97,11 @@ Changes affecting specific commands:
 
     - Include sample name in the output header with `-H` whenever it makes sense (#1992)
 
+* bcftools view
+
+    - Add new `-A, --trim-unseen-allele` option to remove the unseen allele <*> or <NON_REF>
+      at variant sites (`-A`) or all sites (`-AA`)
+
 
 ## Release 1.18 (25th July 2023)
 

bcftools.h

@@ -141,4 +141,16 @@ static inline int bcf_double_test(double d, uint64_t value)
 #define bcf_double_is_missing(x)     bcf_double_test((x),bcf_double_missing)
 #define bcf_double_is_missing_or_vector_end(x)     (bcf_double_test((x),bcf_double_missing) || bcf_double_test((x),bcf_double_vector_end))
 
+static inline int get_unseen_allele(bcf1_t *line)
+{
+    int i;
+    for (i=1; i<line->n_allele; i++)
+    {
+        if ( !strcmp(line->d.allele[i],"<*>") ) return i;
+        if ( !strcmp(line->d.allele[i],"<NON_REF>") ) return i;
+        if ( !strcmp(line->d.allele[i],"<X>") ) return i;
+    }
+    return 0;
+}
+
 #endif

doc/bcftools.1

@@ -2,12 +2,12 @@
 .\"     Title: bcftools
 .\"    Author: [see the "AUTHOR(S)" section]
 .\" Generator: Asciidoctor 2.0.16
-.\"      Date: 2023-10-23
+.\"      Date: 2023-11-09
 .\"    Manual: \ \&
 .\"    Source: \ \&
 .\"  Language: English
 .\"
-.TH "BCFTOOLS" "1" "2023-10-23" "\ \&" "\ \&"
+.TH "BCFTOOLS" "1" "2023-11-09" "\ \&" "\ \&"
 .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\-10\-23 09:31 CEST\fP and refers to bcftools git version \fB1.18\-19\-g6374c6a+\fP.
+This manual page was last updated \fB2023\-11\-09 17:03 GMT\fP and refers to bcftools git version \fB1.18\-24\-g39a81be+\fP.
 .SS "BCF1"
 .sp
 The obsolete BCF1 format output by versions of samtools <= 0.1.19 is \fBnot\fP
@@ -2500,7 +2500,13 @@ in\-memory sorting and DIR is the temporary directory for external sorting. This
 Stop after first record to estimate required time.
 .RE
 .sp
-\fB\-e, \-\-error\-probability\fP \fIINT\fP
+\fB\-e, \-\-exclude\fP [\fIqry\fP|\fIgt\fP]:\*(AqEXPRESSION\*(Aq
+.RS 4
+Exclude sites from query file (\fIqry:\fP) or genotype file (\fIgt:\fP) for which \fIEXPRESSION\fP is true.
+For valid expressions see \fBEXPRESSIONS\fP.
+.RE
+.sp
+\fB\-E, \-\-error\-probability\fP \fIINT\fP
 .RS 4
 Interpret genotypes and genotype likelihoods probabilistically. The value of \fIINT\fP
 represents genotype quality when GT tag is used (e.g. Q=30 represents one error in 1,000 genotypes and
@@ -2510,13 +2516,20 @@ non\-zero integer can be provided).
 .br
 \~
 .br
-If \fB\-e\fP is set to 0, the discordance score can be interpreted as the number of mismatching genotypes,
+If \fB\-E\fP is set to 0, the discordance score can be interpreted as the number of mismatching genotypes,
 but only in the GT\-vs\-GT matching mode. See the \fB\-u, \-\-use\fP option below for additional notes and caveats.
 \~
 .br
 \~
 .br
-If performance is an issue, set \fB\-e 0\fP for faster run times but less accurate results.
+If performance is an issue, set \fB\-E 0\fP for faster run times but less accurate results.
+\~
+.br
+\~
+.br
+Note that in previous versions of bcftools (\(lA1.18), this option used to be a smaller case \fB\-e\fP. It
+changed to make room for the filtering option \fB\-e, \-\-exclude\fP to stay consistent across other
+commands.
 .RE
 .sp
 \fB\-g, \-\-genotypes\fP \fIFILE\fP
@@ -2529,6 +2542,12 @@ VCF/BCF file with reference genotypes to compare against
 Homozygous genotypes only, useful with low coverage data (requires \fB\-g, \-\-genotypes\fP)
 .RE
 .sp
+\fB\-i, \-\-include\fP [\fIqry\fP|\fIgt\fP]:\*(AqEXPRESSION\*(Aq
+.RS 4
+Include sites from query file (\fIqry:\fP) or genotype file (\fIgt:\fP) for which \fIEXPRESSION\fP is true.
+For valid expressions see \fBEXPRESSIONS\fP.
+.RE
+.sp
 \fB\-\-n\-matches\fP \fIINT\fP
 .RS 4
 Print only top INT matches for each sample, 0 for unlimited. Use negative value
@@ -2542,6 +2561,16 @@ Disable calculation of HWE probability to reduce memory requirements with
 comparisons between very large number of sample pairs.
 .RE
 .sp
+\fB\-o, \-\-output\fP \fIFILE\fP
+.RS 4
+Write to \fIFILE\fP rather than to standard output, where it is written by default.
+.RE
+.sp
+\fB\-O, \-\-output\-type\fP \fIt\fP|\fIz\fP
+.RS 4
+Write a plain (\fIt\fP) or compressed (\fIz\fP) text tab\-delimited output.
+.RE
+.sp
 \fB\-p, \-\-pairs\fP \fILIST\fP
 .RS 4
 A comma\-separated list of sample pairs to compare. When the \fB\-g\fP option is given, the first
@@ -2647,6 +2676,14 @@ By default, all header lines are displayed.
 Also display the first \fIINT\fP variant records.
 By default, no variant records are displayed.
 .RE
+.sp
+\fB\-s, \-\-samples\fP \fIINT\fP
+.RS 4
+Display the first \fIINT\fP variant records including the last #CHROM header line with samples.
+Running with \fB\-s 0\fP alone outputs the #CHROM header line only. Note that
+the list of samples, with each sample per line, can be obtained with \f(CRbcftools query\fP using
+the option \fB\-l, \-\-list\-samples\fP.
+.RE
 .SS "bcftools index [\fIOPTIONS\fP] \fIin.bcf\fP|\fIin.vcf.gz\fP"
 .sp
 Creates index for bgzip compressed VCF/BCF files for random access. CSI
@@ -2950,9 +2987,11 @@ maximum number of alternate alleles that can be included in the PL tag. The defa
 is 0 which disables the feature and outputs values for all alternate alleles.
 .RE
 .sp
-\fB\-m, \-\-merge\fP \fIsnps\fP|\fIindels\fP|\fIboth\fP|\fIsnp\-ins\-del\fP|\fIall\fP|\fInone\fP|\fIid\fP
+\fB\-m, \-\-merge\fP \fIsnps\fP|\fIindels\fP|\fIboth\fP|\fIsnp\-ins\-del\fP|\fIall\fP|\fInone\fP|\fIid\fP[,\fI*\fP]
 .RS 4
-The option controls what types of multiallelic records can be created:
+The option controls what types of multiallelic records can be created. If single asterisk
+\fI\fB\fP is appended, the unobserved allele \fI<\fP>\fP or \fI<NON_REF>\fP will be removed at variant sites;
+if two asterisks \fI**\fP are appended, the unobserved allele will be removed all sites.
 .RE
 .sp
 .if n .RS 4
@@ -2962,6 +3001,8 @@ The option controls what types of multiallelic records can be created:
 \-m snps        ..  allow multiallelic SNP records
 \-m indels      ..  allow multiallelic indel records
 \-m both        ..  both SNP and indel records can be multiallelic
+\-m both,*      ..  same as above but remove <*> (or <NON_REF>) from variant sites
+\-m both,**     ..  same as above but remove <*> (or <NON_REF>) at all sites
 \-m all         ..  SNP records can be merged with indel records
 \-m snp\-ins\-del ..  allow multiallelic SNVs, insertions, deletions, but don\*(Aqt mix them
 \-m id          ..  merge by ID
@@ -5144,6 +5185,12 @@ Automatically index the output file
 .RE
 .SS "Subset options:"
 .sp
+\fB\-A, \-\-trim\-unseen\-alleles\fP
+.RS 4
+remove the unseen allele \fI<*>\fP or \fI<NON_REF>\fP at variant sites when the option is given once (\-A) or
+at all sites when the options is given twice (\fI\-AA\fP).
+.RE
+.sp
 \fB\-a, \-\-trim\-alt\-alleles\fP
 .RS 4
 remove alleles not seen in the genotype fields from the ALT column. Note that if no alternate allele

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-10-23 09:31 CEST</strong> and refers to bcftools git version <strong>1.18-19-g6374c6a+</strong>.</p>
+<p>This manual page was last updated <strong>2023-11-09 17:03 GMT</strong> and refers to bcftools git version <strong>1.18-24-g39a81be+</strong>.</p>
 </div>
 </div>
 <div class="sect2">
@@ -2200,19 +2200,29 @@ <h3 id="gtcheck">bcftools gtcheck [<em>OPTIONS</em>] [<strong>-g</strong> <em>ge
 <dd>
 <p>Stop after first record to estimate required time.</p>
 </dd>
-<dt class="hdlist1"><strong>-e, --error-probability</strong> <em>INT</em></dt>
+<dt class="hdlist1"><strong>-e, --exclude</strong> [<em>qry</em>|<em>gt</em>]:'EXPRESSION'</dt>
+<dd>
+<p>Exclude sites from query file (<em>qry:</em>) or genotype file (<em>gt:</em>) for which <em>EXPRESSION</em> is true.
+For valid expressions see <strong><a href="#expressions">EXPRESSIONS</a></strong>.</p>
+</dd>
+<dt class="hdlist1"><strong>-E, --error-probability</strong> <em>INT</em></dt>
 <dd>
 <p>Interpret genotypes and genotype likelihoods probabilistically. The value of <em>INT</em>
 represents genotype quality when GT tag is used (e.g. Q=30 represents one error in 1,000 genotypes and
 Q=40 one error in 10,000 genotypes) and is ignored when PL tag is used (in that case an arbitrary
 non-zero integer can be provided).
 &#160;<br>
 &#160;<br>
-If <strong>-e</strong> is set to 0, the discordance score can be interpreted as the number of mismatching genotypes,
+If <strong>-E</strong> is set to 0, the discordance score can be interpreted as the number of mismatching genotypes,
 but only in the GT-vs-GT matching mode. See the <strong>-u, --use</strong> option below for additional notes and caveats.
 &#160;<br>
 &#160;<br>
-If performance is an issue, set <strong>-e 0</strong> for faster run times but less accurate results.</p>
+If performance is an issue, set <strong>-E 0</strong> for faster run times but less accurate results.
+&#160;<br>
+&#160;<br>
+Note that in previous versions of bcftools (&#8656;1.18), this option used to be a smaller case <strong>-e</strong>. It
+changed to make room for the filtering option <strong>-e, --exclude</strong> to stay consistent across other
+commands.</p>
 </dd>
 <dt class="hdlist1"><strong>-g, --genotypes</strong> <em>FILE</em></dt>
 <dd>
@@ -2222,6 +2232,11 @@ <h3 id="gtcheck">bcftools gtcheck [<em>OPTIONS</em>] [<strong>-g</strong> <em>ge
 <dd>
 <p>Homozygous genotypes only, useful with low coverage data (requires <strong>-g, --genotypes</strong>)</p>
 </dd>
+<dt class="hdlist1"><strong>-i, --include</strong> [<em>qry</em>|<em>gt</em>]:'EXPRESSION'</dt>
+<dd>
+<p>Include sites from query file (<em>qry:</em>) or genotype file (<em>gt:</em>) for which <em>EXPRESSION</em> is true.
+For valid expressions see <strong><a href="#expressions">EXPRESSIONS</a></strong>.</p>
+</dd>
 <dt class="hdlist1"><strong>--n-matches</strong> <em>INT</em></dt>
 <dd>
 <p>Print only top INT matches for each sample, 0 for unlimited. Use negative value
@@ -2233,6 +2248,14 @@ <h3 id="gtcheck">bcftools gtcheck [<em>OPTIONS</em>] [<strong>-g</strong> <em>ge
 <p>Disable calculation of HWE probability to reduce memory requirements with
 comparisons between very large number of sample pairs.</p>
 </dd>
+<dt class="hdlist1"><strong>-o, --output</strong> <em>FILE</em></dt>
+<dd>
+<p>Write to <em>FILE</em> rather than to standard output, where it is written by default.</p>
+</dd>
+<dt class="hdlist1"><strong>-O, --output-type</strong> <em>t</em>|<em>z</em></dt>
+<dd>
+<p>Write a plain (<em>t</em>) or compressed (<em>z</em>) text tab-delimited output.</p>
+</dd>
 <dt class="hdlist1"><strong>-p, --pairs</strong> <em>LIST</em></dt>
 <dd>
 <p>A comma-separated list of sample pairs to compare. When the <strong>-g</strong> option is given, the first
@@ -2339,6 +2362,13 @@ <h4 id="_options">Options:</h4>
 <p>Also display the first <em>INT</em> variant records.
 By default, no variant records are displayed.</p>
 </dd>
+<dt class="hdlist1"><strong>-s, --samples</strong> <em>INT</em></dt>
+<dd>
+<p>Display the first <em>INT</em> variant records including the last #CHROM header line with samples.
+Running with <strong>-s 0</strong> alone outputs the #CHROM header line only. Note that
+the list of samples, with each sample per line, can be obtained with <code>bcftools query</code> using
+the option <strong>-l, --list-samples</strong>.</p>
+</dd>
 </dl>
 </div>
 </div>
@@ -2629,9 +2659,11 @@ <h3 id="merge">bcftools merge [<em>OPTIONS</em>] <em>A.vcf.gz</em> <em>B.vcf.gz<
 maximum number of alternate alleles that can be included in the PL tag. The default value
 is 0 which disables the feature and outputs values for all alternate alleles.</p>
 </dd>
-<dt class="hdlist1"><strong>-m, --merge</strong> <em>snps</em>|<em>indels</em>|<em>both</em>|<em>snp-ins-del</em>|<em>all</em>|<em>none</em>|<em>id</em></dt>
+<dt class="hdlist1"><strong>-m, --merge</strong> <em>snps</em>|<em>indels</em>|<em>both</em>|<em>snp-ins-del</em>|<em>all</em>|<em>none</em>|<em>id</em>[,<em>*</em>]</dt>
 <dd>
-<p>The option controls what types of multiallelic records can be created:</p>
+<p>The option controls what types of multiallelic records can be created. If single asterisk
+<em><strong></em> is appended, the unobserved allele <em>&lt;</strong>&gt;</em> or <em>&lt;NON_REF&gt;</em> will be removed at variant sites;
+if two asterisks <em>**</em> are appended, the unobserved allele will be removed all sites.</p>
 </dd>
 </dl>
 </div>
@@ -2641,6 +2673,8 @@ <h3 id="merge">bcftools merge [<em>OPTIONS</em>] <em>A.vcf.gz</em> <em>B.vcf.gz<
 -m snps        ..  allow multiallelic SNP records
 -m indels      ..  allow multiallelic indel records
 -m both        ..  both SNP and indel records can be multiallelic
+-m both,*      ..  same as above but remove &lt;*&gt; (or &lt;NON_REF&gt;) from variant sites
+-m both,**     ..  same as above but remove &lt;*&gt; (or &lt;NON_REF&gt;) at all sites
 -m all         ..  SNP records can be merged with indel records
 -m snp-ins-del ..  allow multiallelic SNVs, insertions, deletions, but don't mix them
 -m id          ..  merge by ID</pre>
@@ -4506,6 +4540,11 @@ <h4 id="_output_options_2">Output options</h4>
 <h4 id="_subset_options">Subset options:</h4>
 <div class="dlist">
 <dl>
+<dt class="hdlist1"><strong>-A, --trim-unseen-alleles</strong></dt>
+<dd>
+<p>remove the unseen allele <em>&lt;*&gt;</em> or <em>&lt;NON_REF&gt;</em> at variant sites when the option is given once (-A) or
+at all sites when the options is given twice (<em>-AA</em>).</p>
+</dd>
 <dt class="hdlist1"><strong>-a, --trim-alt-alleles</strong></dt>
 <dd>
 <p>remove alleles not seen in the genotype fields from the ALT column. Note that if no alternate allele
@@ -5327,7 +5366,7 @@ <h2 id="_copying">COPYING</h2>
 </div>
 <div id="footer">
 <div id="footer-text">
-Last updated 2023-10-23 09:31:31 +0200
+Last updated 2023-11-09 16:40:20 UTC
 </div>
 </div>
 </body>

doc/bcftools.txt

@@ -2013,13 +2013,17 @@ For "vertical" merge take a look at *<<concat,bcftools concat>>* or *<<norm,bcft
     maximum number of alternate alleles that can be included in the PL tag. The default value
     is 0 which disables the feature and outputs values for all alternate alleles.
 
-*-m, --merge* 'snps'|'indels'|'both'|'snp-ins-del'|'all'|'none'|'id'::
-    The option controls what types of multiallelic records can be created:
+*-m, --merge* 'snps'|'indels'|'both'|'snp-ins-del'|'all'|'none'|'id'[,'*']::
+    The option controls what types of multiallelic records can be created. If single asterisk
+    '*' is appended, the unobserved allele '<*>' or '<NON_REF>' will be removed at variant sites;
+    if two asterisks '**' are appended, the unobserved allele will be removed all sites.
 ----
 -m none        ..  no new multiallelics, output multiple records instead
 -m snps        ..  allow multiallelic SNP records
 -m indels      ..  allow multiallelic indel records
 -m both        ..  both SNP and indel records can be multiallelic
+-m both,*      ..  same as above but remove <*> (or <NON_REF>) from variant sites
+-m both,**     ..  same as above but remove <*> (or <NON_REF>) at all sites
 -m all         ..  SNP records can be merged with indel records
 -m snp-ins-del ..  allow multiallelic SNVs, insertions, deletions, but don't mix them
 -m id          ..  merge by ID
@@ -3397,6 +3401,10 @@ Convert between VCF and BCF. Former *bcftools subset*.
 
 
 ==== Subset options:
+*-A, --trim-unseen-alleles*::
+    remove the unseen allele '<*>' or '<NON_REF>' at variant sites when the option is given once (-A) or
+    at all sites when the options is given twice ('-AA').
+
 *-a, --trim-alt-alleles*::
     remove alleles not seen in the genotype fields from the ALT column. Note that if no alternate allele
     remains after trimming, the record itself is not removed but ALT is set to ".".