Merge branch 'develop' of pd3-github:samtools/bcftools into develop

Author: pd3

NEWS

@@ -6,6 +6,11 @@ Changes affecting specific commands:
 
     - New `-s, --samples` to print per-sample HWE probability, geometric mean
 
+* bcftools conveert
+
+    - Make --gvcf2vcf work around malformed gVCF records where INFO/END overlaps
+      the next record
+
 * bcftools mpileup
 
     - Remove unused experimental INFO/MIN_PL_SUM annotation

test/convert.gvcf.2.out

@@ -0,0 +1,20 @@
+##fileformat=VCFv4.2
+##FILTER=<ID=PASS,Description="All filters passed">
+##reference=file://some/path/human_g1k_v37.fasta
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the region described in this record">
+##contig=<ID=22,length=450>
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	smpl
+22	10	.	C	.	.	.	.	GT	0/0
+22	11	.	T	.	.	.	.	GT	0/0
+22	12	.	T	.	.	.	.	GT	0/0
+22	13	.	G	.	.	.	.	GT	0/0
+22	14	.	G	.	.	.	.	GT	0/0
+22	15	.	C	.	.	.	.	GT	0/0
+22	16	.	C	.	.	.	.	GT	0/0
+22	17	.	A	.	.	.	.	GT	0/0
+22	18	.	A	.	.	.	.	GT	0/0
+22	19	.	G	.	.	.	.	GT	0/0
+22	20	.	T	.	.	.	.	GT	0/0
+22	21	.	C	T	.	.	END=26	GT	0/0
+22	27	.	C	T	.	.	END=42	GT	0/0

test/convert.gvcf.2.vcf

@@ -0,0 +1,9 @@
+##fileformat=VCFv4.2
+##reference=file://some/path/human_g1k_v37.fasta
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the region described in this record">
+##contig=<ID=22,length=450>
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	smpl
+22	10	.	C	.	.	.	END=40	GT	0/0
+22	21	.	C	T	.	.	END=26	GT	0/0
+22	27	.	C	T	.	.	END=42	GT	0/0

test/convert.gvcf.out

@@ -3,7 +3,7 @@
 ##reference=file://some/path/human_g1k_v37.fasta
 ##FORMAT=<ID=GQX,Number=1,Type=Integer,Description="Minimum of {Genotype quality assuming variant position,Genotype quality assuming non-variant position}">
 ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
-##FORMAT=<ID=GQ,Number=1,Type=Float,Description="Genotype Quality">
+##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
 ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Filtered basecall depth used for site genotyping">
 ##FORMAT=<ID=DPF,Number=1,Type=Integer,Description="Basecalls filtered from input prior to site genotyping">
 ##FORMAT=<ID=AD,Number=.,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed. For indels this value only includes reads which confidently support each allele (posterior prob 0.999 or higher that read contains indicated allele vs all other intersecting indel alleles)">
@@ -57,7 +57,7 @@
 ##FORMAT=<ID=OPL,Number=.,Type=Integer,Description="Original PL value before ploidy correction">
 ##INFO=<ID=phastCons,Number=0,Type=Flag,Description="overlaps a phastCons element">
 ##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele, ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/pilot_data/technical/reference/ancestral_alignments/README">
-##INFO=<ID=AF,Number=A,Type=String,Description="1000 Genomes Allele Frequency based on AC/AN; Format: Allele:AlleleFrequency">
+##INFO=<ID=AF,Number=A,Type=Float,Description="1000 Genomes Allele Frequency based on AC/AN; Format: Allele:AlleleFrequency">
 ##INFO=<ID=AMR_AF,Number=A,Type=String,Description="1000 Genomes Allele Frequency for samples from AMR population based on AC/AN; Format: Allele:AlleleFrequency">
 ##INFO=<ID=ASN_AF,Number=A,Type=String,Description="1000 Genomes Allele Frequency for samples from ASN population based on AC/AN; Format: Allele:AlleleFrequency">
 ##INFO=<ID=AFR_AF,Number=A,Type=String,Description="1000 Genomes Allele Frequency for samples from AFR population based on AC/AN; Format: Allele:AlleleFrequency">

test/convert.gvcf.vcf

@@ -3,7 +3,7 @@
 ##reference=file://some/path/human_g1k_v37.fasta
 ##FORMAT=<ID=GQX,Number=1,Type=Integer,Description="Minimum of {Genotype quality assuming variant position,Genotype quality assuming non-variant position}">
 ##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
-##FORMAT=<ID=GQ,Number=1,Type=Float,Description="Genotype Quality">
+##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
 ##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Filtered basecall depth used for site genotyping">
 ##FORMAT=<ID=DPF,Number=1,Type=Integer,Description="Basecalls filtered from input prior to site genotyping">
 ##FORMAT=<ID=AD,Number=.,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed. For indels this value only includes reads which confidently support each allele (posterior prob 0.999 or higher that read contains indicated allele vs all other intersecting indel alleles)">
@@ -57,7 +57,7 @@
 ##FORMAT=<ID=OPL,Number=.,Type=Integer,Description="Original PL value before ploidy correction">
 ##INFO=<ID=phastCons,Number=0,Type=Flag,Description="overlaps a phastCons element">
 ##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele, ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/pilot_data/technical/reference/ancestral_alignments/README">
-##INFO=<ID=AF,Number=A,Type=String,Description="1000 Genomes Allele Frequency based on AC/AN; Format: Allele:AlleleFrequency">
+##INFO=<ID=AF,Number=A,Type=Float,Description="1000 Genomes Allele Frequency based on AC/AN; Format: Allele:AlleleFrequency">
 ##INFO=<ID=AMR_AF,Number=A,Type=String,Description="1000 Genomes Allele Frequency for samples from AMR population based on AC/AN; Format: Allele:AlleleFrequency">
 ##INFO=<ID=ASN_AF,Number=A,Type=String,Description="1000 Genomes Allele Frequency for samples from ASN population based on AC/AN; Format: Allele:AlleleFrequency">
 ##INFO=<ID=AFR_AF,Number=A,Type=String,Description="1000 Genomes Allele Frequency for samples from AFR population based on AC/AN; Format: Allele:AlleleFrequency">

test/test.pl

@@ -949,6 +949,7 @@
 run_test(\&test_vcf_convert_hs2vcf,$opts,h=>'convert.hs.gt.hap',s=>'convert.hs.gt.samples',out=>'convert.gt.noHead.vcf',args=>'--hapsample2vcf');
 run_test(\&test_vcf_convert_hs2vcf,$opts,h=>'convert.hs.gt.ids.hap',s=>'convert.hs.gt.samples',out=>'convert.gt.noHead.ids.vcf',args=>'--vcf-ids --hapsample2vcf');
 run_test(\&test_vcf_convert_gvcf,$opts,in=>'convert.gvcf',out=>'convert.gvcf.out',fa=>'gvcf.fa',args=>'--gvcf2vcf -i\'FILTER="PASS"\'');
+run_test(\&test_vcf_convert_gvcf,$opts,in=>'convert.gvcf.2',out=>'convert.gvcf.2.out',fa=>'gvcf.fa',args=>'--gvcf2vcf');
 run_test(\&test_vcf_convert_tsv2vcf,$opts,in=>'convert.23andme',out=>'convert.23andme.vcf',args=>'-c ID,CHROM,POS,AA -s SAMPLE1',fai=>'23andme');
 run_test(\&test_vcf_convert_tsv2vcf,$opts,in=>'convert.tsv',out=>'convert.tsv.vcf',args=>'-c -,CHROM,POS,REF,ALT',fai=>'23andme');
 run_test(\&test_vcf_consensus,$opts,in=>'consensus.gvcf-missing.1',out=>'consensus.gvcf-missing.1.out',fa=>'consensus.gvcf-missing.fa',args=>'--missing N');

vcfconvert.c

@@ -1,6 +1,6 @@
 /*  vcfconvert.c -- convert between VCF/BCF and related formats.
 
-    Copyright (C) 2013-2025 Genome Research Ltd.
+    Copyright (C) 2013-2026 Genome Research Ltd.
 
     Author: Petr Danecek <pd3@sanger.ac.uk>
 
@@ -44,6 +44,7 @@ THE SOFTWARE.  */
 #include "filter.h"
 #include "convert.h"
 #include "tsv2vcf.h"
+#include "vcfbuf.h"
 
 // Logic of the filters: include or exclude sites which match the filters?
 #define FLT_INCLUDE 1
@@ -52,6 +53,7 @@ THE SOFTWARE.  */
 typedef struct _args_t args_t;
 struct _args_t
 {
+    vcfbuf_t *vcfbuf;
     faidx_t *ref;
     filter_t *filter;
     char *filter_str;
@@ -64,7 +66,8 @@ struct _args_t
         int total, skipped, hom_rr, het_ra, hom_aa, het_aa, missing, written;
     } n;
     kstring_t str;
-    int32_t *gts;
+    int32_t *gts, *itmp;
+    int nitmp;
     float *flt;
     int rev_als, output_vcf_ids, hap2dip, gen_3N6;
     int nsamples, *samples, sample_is_file, targets_is_file, regions_is_file, output_type;
@@ -1512,6 +1515,50 @@ static void vcf_to_vcf(args_t *args)
     if ( hts_close(out_fh)!=0 ) error("[%s] Error: close failed .. %s\n", __func__,args->outfname);
 }
 
+static bcf1_t *gvcf_next_line(args_t *args, int *end1)
+{
+    *end1 = 0;
+    while ( vcfbuf_nsites(args->vcfbuf)!=2 && bcf_sr_next_line(args->files) )
+    {
+        bcf1_t *rec = bcf_sr_get_line(args->files,0);
+        args->files->readers[0].buffer[0] = vcfbuf_push(args->vcfbuf, rec);
+    }
+    bcf1_t *rec = vcfbuf_flush(args->vcfbuf, 0);
+    if ( !rec ) return NULL;
+
+    // is it a gVCF record?
+    // - ALT must be one of ., <*>, <X>, <NON_REF>
+    // - INFO/END must be present
+    int i, gallele = -1;
+    if ( rec->n_allele==1 ) gallele = 0; // illumina/bcftools-call gvcf (if INFO/END present)
+    else if ( rec->d.allele[1][0]=='<' )
+    {
+        for (i=1; i<rec->n_allele; i++)
+        {
+            if ( rec->d.allele[i][1]=='*' && rec->d.allele[i][2]=='>' && rec->d.allele[i][3]=='\0' ) { gallele = i; break; } // mpileup/spec compliant gVCF
+            if ( rec->d.allele[i][1]=='X' && rec->d.allele[i][2]=='>' && rec->d.allele[i][3]=='\0' ) { gallele = i; break; } // old mpileup gVCF
+            if ( strcmp(rec->d.allele[i],"<NON_REF>")==0 ) { gallele = i; break; }               // GATK gVCF
+        }
+    }
+    if ( gallele<0 ) return rec;
+
+    int nend = bcf_get_info_int32(args->header,rec,"END",&args->itmp,&args->nitmp);
+    if ( nend!=1 ) return rec;
+    *end1 = args->itmp[0];
+
+    bcf1_t *peek = vcfbuf_peek(args->vcfbuf, 0);
+    if ( peek && rec->rid==peek->rid && (*end1)-1>=peek->pos )
+    {
+        static int warned = 0;
+        if ( !warned )
+        {
+            hts_log_warning("Malformed gVCF: INFO/END at %s:%"PRIhts_pos" overlaps the next record",bcf_seqname(args->header,rec),rec->pos+1);
+            warned = 1;
+        }
+        *end1 = rec->pos < peek->pos ? peek->pos : 0;
+    }
+    return rec;
+}
 static void gvcf_to_vcf(args_t *args)
 {
     if ( !args->ref_fname ) error("--gvcf2vcf requires the --fasta-ref option\n");
@@ -1529,15 +1576,16 @@ static void gvcf_to_vcf(args_t *args)
     bcf_hdr_t *hdr = bcf_sr_get_header(args->files,0);
     if (args->record_cmd_line) bcf_hdr_append_version(hdr, args->argc, args->argv, "bcftools_convert");
     if ( bcf_hdr_write(out_fh,hdr)!=0 ) error("[%s] Error: cannot write to %s\n", __func__,args->outfname);
-    if ( init_index2(out_fh,hdr,args->outfname,&args->index_fn,
-                     args->write_index)<0 )
+    if ( init_index2(out_fh,hdr,args->outfname,&args->index_fn,args->write_index)<0 )
         error("Error: failed to initialise index for %s\n",args->outfname);
 
-    int32_t *itmp = NULL, nitmp = 0;
+    args->vcfbuf = vcfbuf_init(hdr, 0);
+    vcfbuf_set(args->vcfbuf,VCFBUF_DUMMY,1);
 
-    while ( bcf_sr_next_line(args->files) )
+    int end1 = 0;
+    bcf1_t *line;
+    while ( (line=gvcf_next_line(args, &end1)) )
     {
-        bcf1_t *line = bcf_sr_get_line(args->files,0);
         if ( args->filter )
         {
             int pass = filter_test(args->filter, line, NULL);
@@ -1549,39 +1597,15 @@ static void gvcf_to_vcf(args_t *args)
             }
         }
 
-        // check if alleles compatible with being a gVCF record
-        // ALT must be one of ., <*>, <X>, <NON_REF>
-        // check for INFO/END is below
-        int i, gallele = -1;
-        if (line->n_allele==1)
-            gallele = 0; // illumina/bcftools-call gvcf (if INFO/END present)
-        else if ( line->d.allele[1][0]=='<' )
-        {
-            for (i=1; i<line->n_allele; i++)
-            {
-                if ( line->d.allele[i][1]=='*' && line->d.allele[i][2]=='>' && line->d.allele[i][3]=='\0' ) { gallele = i; break; } // mpileup/spec compliant gVCF
-                if ( line->d.allele[i][1]=='X' && line->d.allele[i][2]=='>' && line->d.allele[i][3]=='\0' ) { gallele = i; break; } // old mpileup gVCF
-                if ( strcmp(line->d.allele[i],"<NON_REF>")==0 ) { gallele = i; break; }               // GATK gVCF
-            }
-        }
-
-        // no gVCF compatible alleles
-        if (gallele<0)
+        if ( !end1 )
         {
             if ( bcf_write(out_fh,hdr,line)!=0 ) error("[%s] Error: cannot write to %s\n", __func__,args->outfname);
             continue;
         }
 
-        int nend = bcf_get_info_int32(hdr,line,"END",&itmp,&nitmp);
-        if ( nend!=1 )
-        {
-            // No INFO/END => not gVCF record
-            if ( bcf_write(out_fh,hdr,line)!=0 ) error("[%s] Error: cannot write to %s\n", __func__,args->outfname);
-            continue;
-        }
         bcf_update_info_int32(hdr,line,"END",NULL,0);
         int pos, len;
-        for (pos=line->pos; pos<itmp[0]; pos++)
+        for (pos=line->pos; pos<end1; pos++)
         {
             line->pos = pos;
             char *ref = faidx_fetch_seq(args->ref, (char*)bcf_hdr_id2name(hdr,line->rid), line->pos, line->pos, &len);
@@ -1593,7 +1617,7 @@ static void gvcf_to_vcf(args_t *args)
         }
     }
     if ( args->files->errnum ) error("Error: %s\n", bcf_sr_strerror(args->files->errnum));
-    free(itmp);
+    free(args->itmp);
     if ( args->write_index )
     {
         if ( bcf_idx_save(out_fh)<0 )
@@ -1604,6 +1628,7 @@ static void gvcf_to_vcf(args_t *args)
         free(args->index_fn);
     }
     if ( hts_close(out_fh)!=0 ) error("[%s] Error: close failed .. %s\n", __func__,args->outfname);
+    vcfbuf_destroy(args->vcfbuf);
 }
 
 static void usage(void)