Fix a bug for incorrectly formatted gVCF files

When the gVCF contains overlapping blocks, they would trigger
an infinite loop in the program and it would never finish

Resolves #2410

Author: pd3

NEWS

@@ -67,6 +67,9 @@ Changes affecting specific commands:
     - Do not merge symbolic alleles unless they have not just the same type, eg. <DEL>,
       but also length, i.e the INFO/END coordinate (#2362)
 
+    - Fix a bug where an incorrectly formatted gVCF file with overlapping blocks would trigger
+      an infinite loop in the program (#2410)
+
 * bcftools mpileup
 
     - The -r/-R option newly merge overlapping regions, preventing the output of duplicate sites

test/merge.broken-gvcf.1.out

@@ -0,0 +1,9 @@
+##fileformat=VCFv4.2
+##FILTER=<ID=PASS,Description="All filters passed">
+##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	SB
+chr1	80000	id1	G	.	0	PASS	END=90000	GT:CN	0/0:2	0/0:2
+chr1	80000	id2	G	.	15.43	PASS	.	GT:CN	0/0:2	0/0:2

test/merge.broken-gvcf.a.vcf

@@ -0,0 +1,8 @@
+##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

test/merge.broken-gvcf.b.vcf

@@ -0,0 +1,8 @@
+##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

test/test.pl

@@ -61,6 +61,7 @@
 run_test(\&test_vcf_isec,$opts,in=>['isec-miss.1.1','isec-miss.1.2','isec-miss.1.3'],out=>'isec-miss.1.1.out',args=>'-R {PATH}/isec-miss.1.regs.txt -n +1');
 run_test(\&test_vcf_isec,$opts,in=>['isec-miss.2.1','isec-miss.2.2','isec-miss.2.3'],out=>'isec-miss.2.1.out',args=>'-n +1 -r 20:100,20:140,12:55,20:140,20:100');
 run_test(\&test_vcf_isec,$opts,in=>['isec-miss.2.1','isec-miss.2.2','isec-miss.2.3'],out=>'isec-miss.2.1.out',args=>'-R {PATH}/isec-miss.1.regs.txt -n +1');
+run_test(\&test_vcf_merge,$opts,in=>['merge.broken-gvcf.a','merge.broken-gvcf.b'],out=>'merge.broken-gvcf.1.out',args=>'');
 run_test(\&test_vcf_merge,$opts,in=>['merge.symbolic.1.a','merge.symbolic.1.b'],out=>'merge.symbolic.1.1.out',args=>'');
 run_test(\&test_vcf_merge,$opts,in=>['merge.multiallelics.1.a','merge.multiallelics.1.b'],out=>'merge.multiallelics.1.1.out',args=>'--merge none');
 run_test(\&test_vcf_merge,$opts,in=>['merge.multiallelics.1.a','merge.multiallelics.1.b'],out=>'merge.multiallelics.1.1.out',args=>'--merge both');

vcfmerge.c

@@ -2723,20 +2723,31 @@ void gvcf_flush(args_t *args, int done)
     }
 }
 
-static inline int is_gvcf_block(bcf1_t *line)
+static inline int is_gvcf_block(args_t *args, bcf1_t *line)
 {
     if ( line->rlen<=1 ) return 0;
     if ( strlen(line->d.allele[0])==line->rlen ) return 0;
-    if ( line->n_allele==1 ) return 1;
+    if ( line->n_allele==1 ) goto is_gvcf;
 
     int i;
     for (i=1; i<line->n_allele; i++)
     {
-        if ( !strcmp(line->d.allele[i],"<*>") ) return 1;
-        if ( !strcmp(line->d.allele[i],"<NON_REF>") ) return 1;
-        if ( !strcmp(line->d.allele[i],"<X>") ) return 1;
+        if ( !strcmp(line->d.allele[i],"<*>") ) goto is_gvcf;
+        if ( !strcmp(line->d.allele[i],"<NON_REF>") ) goto is_gvcf;
+        if ( !strcmp(line->d.allele[i],"<X>") ) goto is_gvcf;
     }
     return 0;
+
+is_gvcf:
+    maux_t *ma = args->maux;
+    if ( !ma->gvcf )
+    {
+        args->do_gvcf = 1;
+        ma->gvcf = (gvcf_aux_t*) calloc(ma->n,sizeof(gvcf_aux_t));  // -Walloc-size-larger-than gives a harmless warning caused by signed integer ma->n
+        for (i=0; i<ma->n; i++)
+            ma->gvcf[i].line = bcf_init1();
+    }
+    return 1;
 }
 
 /*
@@ -2776,7 +2787,7 @@ void gvcf_stage(args_t *args, int pos)
         int irec = maux->buf[i].beg;
         bcf_hdr_t *hdr = bcf_sr_get_header(files, i);
         bcf1_t *line = args->files->readers[i].buffer[irec];
-        int ret = is_gvcf_block(line) ? bcf_get_info_int32(hdr,line,"END",&end,&nend) : 0;
+        int ret = is_gvcf_block(args,line) ? bcf_get_info_int32(hdr,line,"END",&end,&nend) : 0;
         if ( ret==1 )
         {
             if ( end[0] == line->pos + 1 )  // POS and INFO/END are identical, treat as if a normal w/o INFO/END
@@ -3133,7 +3144,7 @@ int can_merge(args_t *args)
                     var_type &= ~VCF_INDEL;
                 }
                 var_type = var_type ? var_type<<1 : ref_mask;
-                if ( args->do_gvcf && is_gvcf_block(line) ) var_type |= ref_mask;
+                if ( args->do_gvcf && is_gvcf_block(args,line) ) var_type |= ref_mask;
                 buf->rec[j].var_types = var_type;
             }
             maux->var_types |= buf->rec[j].var_types;
@@ -3233,7 +3244,8 @@ void stage_line(args_t *args)
             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 && is_gvcf_block(args,buf->lines[j]) ) is_gvcf = 1;
+                if ( is_gvcf && buf->rec[j].skip && !maux->gvcf[i].active ) continue;
                 if ( !is_gvcf ) continue;   // done or not compatible
             }
             if ( args->merge_by_id ) break;     // if merging by ID and the line is compatible, the this is THE line