Fix a bug when update of INFO/END results in assertion error

When INFO/END is updated from a text annotation file, VCF record that overlaps
multiple annotation records could update rlen on the fly and lead to negative
overlap length calculation.

Fixes #1957

Author: pd3

NEWS

@@ -14,6 +14,8 @@ Changes affecting specific commands:
 
     - Fix a bug where the `-m` function did not respect the `--min-overlap` option (#1869)
 
+    - Fix a bug when update of INFO/END results in assertion error (#1957)
+
 * bcftools concat
 
     - New option `--drop-genotypes`

test/annotate34.out

@@ -0,0 +1,14 @@
+##fileformat=VCFv4.2
+##FILTER=<ID=PASS,Description="All filters passed">
+##contig=<ID=chr21,length=45090682>
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
+##FILTER=<ID=HET1,Description="Heterozygous in the first haplotype">
+##FILTER=<ID=HET2,Description="Heterozygous in the second haplotype">
+##INFO=<ID=SVTYPE,Number=1,Type=String,Description="SVTYPE">
+##INFO=<ID=SVLEN,Number=1,Type=Integer,Description="SVLEN">
+##INFO=<ID=END,Number=1,Type=Integer,Description="End coordinate in reference for SV">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	syndip
+chr21	8914240	.	gttccattccattccattcaattccattccattgcattccattccattccattcca	G	30	HET2	SVTYPE=DEL;SVLEN=55;END=8914295	GT:AD	0|.:3,1
+chr21	8914680	.	tattccattccattcc	TATTCCATTCCATTCCATTCCATTCCATTCCATTCCATTCCATTCTAGTTGATTCCATTCCATTCCATCCCGTTCCATTCCATTCCGTTACTTTCTATTCCATTCCATTCCATTCC	30	HET2	SVTYPE=INS;SVLEN=100;END=8914681	GT:AD	0|.:1,1
+chr21	8914690	.	c	CATTCCATTCCATTCCATTCCATTCTAGTTGATTCCATTCCATTCCATCCCGTTCCATTCCATTCCGTTACTTTCT	30	HET2	SVTYPE=INS;SVLEN=75;END=8914691	GT:AD	0|.:2,1

test/annotate34.vcf

@@ -0,0 +1,13 @@
+##fileformat=VCFv4.2
+##FILTER=<ID=PASS,Description="All filters passed">
+##contig=<ID=chr21,length=45090682>
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
+##FILTER=<ID=HET1,Description="Heterozygous in the first haplotype">
+##FILTER=<ID=HET2,Description="Heterozygous in the second haplotype">
+##INFO=<ID=SVTYPE,Number=1,Type=String,Description="SVTYPE">
+##INFO=<ID=SVLEN,Number=1,Type=Integer,Description="SVLEN">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	syndip
+chr21	8914240	.	gttccattccattccattcaattccattccattgcattccattccattccattcca	G	30	HET2	SVTYPE=DEL;SVLEN=55	GT:AD	0|.:3,1
+chr21	8914680	.	tattccattccattcc	TATTCCATTCCATTCCATTCCATTCCATTCCATTCCATTCCATTCTAGTTGATTCCATTCCATTCCATCCCGTTCCATTCCATTCCGTTACTTTCTATTCCATTCCATTCCATTCC	30	HET2	SVTYPE=INS;SVLEN=100	GT:AD	0|.:1,1
+chr21	8914690	.	c	CATTCCATTCCATTCCATTCCATTCTAGTTGATTCCATTCCATTCCATCCCGTTCCATTCCATTCCGTTACTTTCT	30	HET2	SVTYPE=INS;SVLEN=75	GT:AD	0|.:2,1

test/annots34.tab

@@ -0,0 +1,3 @@
+chr21	8914240	8914240	8914295
+chr21	8914680	8914680	8914681
+chr21	8914690	8914690	8914691

test/test.pl

@@ -523,6 +523,7 @@
 run_test(\&test_vcf_annotate,$opts,in=>'annotate28',tab=>'annots28',out=>'annotate28.3.out',args=>'-c CHROM,POS,REF,ALT,FMT/TEST -s smpl1');
 run_test(\&test_vcf_annotate,$opts,in=>'annotate28',tab=>'annots28',out=>'annotate28.4.out',args=>'-c CHROM,POS,REF,ALT,FMT/TEST -s smpl2');
 run_test(\&test_vcf_annotate,$opts,in=>'annotate',out=>'annotate.33.out',args=>'-m XXX');
+run_test(\&test_vcf_annotate,$opts,in=>'annotate34',tab=>'annots34',out=>'annotate34.out',args=>q[-c CHROM,FROM,TO,INFO/END -H '##INFO=<ID=END,Number=1,Type=Integer,Description="End coordinate in reference for SV">']);
 run_test(\&test_vcf_plugin,$opts,in=>'checkploidy',out=>'checkploidy.out',cmd=>'+check-ploidy --no-version');
 run_test(\&test_vcf_plugin,$opts,in=>'checkploidy.2',out=>'checkploidy.2.out',cmd=>'+check-ploidy --no-version');
 run_test(\&test_vcf_plugin,$opts,in=>'checkploidy.2',out=>'checkploidy.3.out',cmd=>'+check-ploidy --no-version',args=>'-- -m');

vcfannotate.c

@@ -3094,6 +3094,7 @@ static void annotate(args_t *args, bcf1_t *line)
         for (j=0; j<args->ncols; j++) args->cols[j].done = 0;
         if ( regidx_overlap(args->tgt_idx, bcf_seqname(args->hdr,line),line->pos,line->pos+line->rlen-1, args->tgt_itr) )
         {
+            hts_pos_t vcf_end = line->pos + line->rlen - 1;
             while ( regitr_overlap(args->tgt_itr) )
             {
                 annot_line_t *tmp = &args->alines[0];
@@ -3104,7 +3105,7 @@ static void annotate(args_t *args, bcf1_t *line)
                 // Check min overlap
                 int len_ann = tmp->end - tmp->start + 1;
                 int len_vcf = line->rlen;
-                int isec = (tmp->end < line->pos+line->rlen-1 ? tmp->end : line->pos+line->rlen-1) - (tmp->start > line->pos ? tmp->start : line->pos) + 1;
+                int isec = (tmp->end < vcf_end ? tmp->end : vcf_end) - (tmp->start > line->pos ? tmp->start : line->pos) + 1;
                 assert( isec > 0 );
                 if ( args->min_overlap_ann && args->min_overlap_ann > (float)isec/len_ann ) continue;
                 if ( args->min_overlap_vcf && args->min_overlap_vcf > (float)isec/len_vcf ) continue;