END field-related bcftools merge infinite loop and infinite output

[jmarshall]

Given the following two VCF files (which will then be bgzipped and tabix-indexed):

##fileformat=VCFv4.2
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=CN,Number=1,Type=Integer,Description="Copy number">
##INFO=<ID=END,Number=1,Type=Integer,Description="End coordinate of the variant">
##contig=<ID=chr1,length=248956422>
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	SA
chr1	80000	id1	G	.	0	PASS	END=90000	GT:CN	0/0:2
chr1	80000	id2	G	.	12.34	PASS	END=95000	GT:CN	0/0:2

and

##fileformat=VCFv4.2
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=CN,Number=1,Type=Integer,Description="Copy number">
##INFO=<ID=END,Number=1,Type=Integer,Description="End coordinate of the variant">
##contig=<ID=chr1,length=248956422>
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	SB
chr1	80000	id1	G	.	0	PASS	END=90000	GT:CN	0/0:2
chr1	80000	id2	G	.	15.43	PASS	END=95000	GT:CN	0/0:2

(attached as infinite.a.vcf.txt and infinite.b.vcf.txt)

the following bcftools merge command runs forever generating the same output line over and over:

bgzip infinite.a.vcf
bgzip infinite.b.vcf
tabix infinite.a.vcf.gz
tabix infinite.b.vcf.gz

bcftools merge infinite.a.vcf.gz infinite.b.vcf.gz -o out.vcf -m all -0

This outputs the merged first record repeatedly forever:

chr1	80000	id1	G	.	0	PASS	END=90000	GT:CN	0/0:2	0/0:2
chr1	80000	id1	G	.	0	PASS	END=90000	GT:CN	0/0:2	0/0:2
chr1	80000	id1	G	.	0	PASS	END=90000	GT:CN	0/0:2	0/0:2
chr1	80000	id1	G	.	0	PASS	END=90000	GT:CN	0/0:2	0/0:2
[…]

BCFtools 1.17 produced appropriately merged output for these input files. Bisecting shows that the issue arrived with 60c65eb.

[jkbonfield]

I took a look and see it's looping in this:

        while ( can_merge(args) )                                                                                        
        {                                                                                                                
            stage_line(args);                                                                                            
            merge_line(args);                                                                                            
        }                                                                                                                

The problem here is in theory as it merges it's meant to strike them out so can_merge ultimately finds nothing left to do, but this isn't happening. I tried to follow the changes in the commit and got lost. There's over 500 lines of code changed (of which >400 are entirely new). This needs @pd3's knowledge.

The good news is if it's taken this long (over 2 years and 4 releases) to be spotted, then I don't see it as a critical regression that holds up the next release if we fail to find a solution quickly.

[jmarshall]

Indeed it is; if you turn debug_maux(args) on, you can see that on the first iteration it marks the first lines as DONE but it never manages to get the second lines marked as DONE:

Alleles to merge at 80000, nals=1
 reader 0 (k=0-2):      !G      G
 reader 1 (k=0-2):      !G      G
 counts:    4x G

Alleles to merge at 80000, nals=1
 reader 0 (k=0-2):  DONE        G
 reader 1 (k=0-2):  DONE        G
 counts:    2x G

Alleles to merge at 80000, nals=1
 reader 0 (k=0-2):  DONE        G
 reader 1 (k=0-2):  DONE        G
 counts:    2x G

[…forever…]

For …reasons…, the pipeline we observed this on has been pinned to bcftools-1.16 (prior to that large commit, which landed in 1.18) and we have only recently tried it with more recent bcftools. Hence why we hadn't spotted this before. Somewhat surprising that others have not also reported this.

[jkbonfield]

Ah thanks for the hint. I wasn't aware of that debugging.

Yes, I can see the culprit in stage_line https://github.com/samtools/bcftools/blob/develop/vcfmerge.c#L3233-L3238

            if ( buf->rec[j].skip )
            {
                int is_gvcf = maux->gvcf && maux->gvcf[i].active ? 1 : 0;
                if ( !is_gvcf && is_gvcf_block(buf->lines[j]) ) is_gvcf = 1;
                if ( !is_gvcf ) continue;   // done or not compatible
            }

So being marked as DONE (skip true) isn't sufficient to skip it and continue to the next entry. In this case as it's marked as gVCF the continue statement doesn't happen. This causes it to repeatedly attempt to merge in the same record with j==0 instead of moving on to the next one.

The line used to just be if ( buf->rec[j].skip ) continue. However this was deliberately changed, obviously, and the commit refers to fixing gVCF merging.

It looks like is_gvcf_block is claiming it's gVCF based purely on if ( line->n_allele==1) return 1. I'm guessing it's thinking about symbolic alleles <*> and <NON_REF> (also listed explicitly), but the spec doesn't say anything about unspecified allelies being equivalent to gVCF also the spec does cryptically mention them in the same section "Representing unspecified alleles and REF-only blocks (gVCF)". Whether that is meant to imply an unspecified allele is gVCF is unknown. Edit: the Illumina documentation on gVCF does use ALT ".", but it's excluded from the formal VCF specification. This is, I assume, an error.

That said, even if that does make it gVCF, it's not clear how this would ever avoid the infinite loop.

[jkbonfield]

Perhaps the issue is that maux->gvcf is only set when using -g and setting args->do_gvcf. Without that, we shouldn't be checking for the active gvcf status.

Also I still don't get it. The order of the lines if maux->gvcf && maux->gvcf[i].active when we're are gvcf, otherwise not. However if we're not gVCF it then reverts back to being one if is_gvcf_block finds <*> etc in the alleles. So even if we're gone through the gVCF blocks and set active to 0 then it'll get reinstated again by the following line. These feel back to front.

That said, I don't really understand what merging gVCF means when we have lots of different sized blocks. Do we end up doing some complex combinatorial intersection and breaking bigger blocks into smaller ones based on whatever sample ends next?

Edit: on your real data, does bcftools merge -g - -m all in1.vcf.gz in2.vcf.gz -o out.vcf.gz not do what you want? Ie explicitly telling it to do the gVCF matching? It does rather feel like this is a bug caused by assuming gVCF mode is enabled during merges of gVCF like records.

[jmarshall]

Yes, the files this is distilled from are actually gVCFs really. I hadn't spotted that about the record data triggering a gVCF mode in the code.

I tried it with -g GRCh38.fa too, and the behaviour was better (no infinite loop) but IIRC the output was differently incorrect.

[pd3]

Thank you for the issue. The program seems to be confused by the conflicting information given in the VCF - there are two overlapping gVCF blocks, which should never happen. A quick workaround is to run bcftools norm -d exact. I'll take a look to see how to resolve situations like this

[jkbonfield]

On Wed, May 28, 2025 at 02:08:25AM -0700, Petr Danecek wrote: Thank you for the issue. The program seems to be confused by the conflicting information given in the VCF - there are two overlapping gVCF blocks, which should never happen.

Thanks Petr. As an aside: it would be useful if statements like this were in the VCF specification. It's extremely vague (pretty much absent!) as to what actually constitutes a gVCF.

…

---------------------------------------------------------------------- The Wellcome Sanger Institute is operated by Genome Research Limited, a charity registered in England with number 1021457 and a company registered in England with number 2742969, whose registered office is Wellcome Sanger Institute, Wellcome Genome Campus, Hinxton, CB10 1SA.

[pd3]

I believe this is now fixed. Thank you for the tests and analysis

[jmarshall]

Spectacular — Thank you, Petr!

The original files are from a little-used part of GATK, and yes the overlapping gVCF blocks are “interesting”. In most cases in the original the overlapping blocks have different FILTER fields, which distinguishes them somewhat but it's still pretty horrible data.

Using -g - made these test files terminate properly, but the full production files the tests were distilled from still caused an infinite loop a few records later.

With @pd3's fix, merging the full production files terminates properly both with and without -g -. Without -g - it produces somewhat different output from bcftools 1.17. However with -g - (which probably should be being used with this gVCF data), bcftools with this fix produces identical output to bcftools 1.17 on these production files. 🎉