Don't expand REF when symbolic <DEL> allele is present

Symbolic <DEL> alleles caused norm to expand REF to the full length of the deletion.
This was not intended and was problematic for long deletions, the REF allele should list
one base only.

Resolves #2029

Author: pd3

NEWS

@@ -40,6 +40,10 @@ Changes affecting specific commands:
 
     - Allow combining -m and -a with --old-rec-tag (#2020)
 
+    - Symbolic <DEL> alleles caused norm to expand REF to the full length of the deletion.
+      This was not intended and problematic for long deletions, the REF allele should list
+      one base only (#2029)
+
 * bcftools query
 
     - Add new `-N, --disable-automatic-newline` option for pre-1.18 query formatting behavior

test/norm.symbolic.1.out

@@ -6,6 +6,6 @@
 ##INFO=<ID=ORI,Number=1,Type=String,Description="Original variant. Format: CHR|POS|REF|ALT|USED_ALT_IDX">
 #CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
 20	15	.	TAC	T	.	.	ORI=20|24|ACA|A
-20	15	.	TAC	<DEL>	.	.	END=17;SVTYPE=DEL;ORI=20|24|A|<DEL>
+20	15	.	T	<DEL>	.	.	END=17;SVTYPE=DEL;ORI=20|24|A|<DEL>
 20	93	.	CAAA	C	.	.	ORI=20|98|AAAA|A
-20	93	.	CAAA	<DEL>	.	.	END=96;SVTYPE=DEL;ORI=20|98|A|<DEL>
+20	93	.	C	<DEL>	.	.	END=96;SVTYPE=DEL;ORI=20|98|A|<DEL>

test/norm.symbolic.2.out

@@ -0,0 +1,10 @@
+##fileformat=VCFv4.2
+##FILTER=<ID=PASS,Description="All filters passed">
+##contig=<ID=20,length=2147483647>
+##INFO=<ID=END,Number=1,Type=Integer,Description="End position">
+##INFO=<ID=SVTYPE,Number=1,Type=String,Description="Type of structural variant">
+##INFO=<ID=ORI,Number=1,Type=String,Description="Original variant. Format: CHR|POS|REF|ALT|USED_ALT_IDX">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
+20	15	.	TAC	T,<DEL>	.	.	END=17;SVTYPE=DEL;ORI=20|24|ACA|A,<DEL>
+20	93	.	CAAA	C	.	.	ORI=20|98|AAAA|A
+20	93	.	C	<DEL>	.	.	END=96;SVTYPE=DEL;ORI=20|98|A|<DEL>

test/norm.symbolic.2.vcf

@@ -0,0 +1,8 @@
+##fileformat=VCFv4.2
+##contig=<ID=20,length=2147483647>
+##INFO=<ID=END,Number=1,Type=Integer,Description="End position">
+##INFO=<ID=SVTYPE,Number=1,Type=String,Description="Type of structural variant">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
+20	24	.	ACA	A,<DEL>	.	.	END=26;SVTYPE=DEL
+20	98	.	AAAA	A	.	.	.
+20	98	.	A	<DEL>	.	.	END=101;SVTYPE=DEL

test/test.pl

@@ -289,6 +289,7 @@
 run_test(\&test_vcf_norm,$opts,in=>'norm.phased-split',out=>'norm.phased-split.1.out',args=>'-m -any');
 run_test(\&test_vcf_norm,$opts,in=>'norm.phased-join',out=>'norm.phased-join.1.out',args=>'-m +any');
 run_test(\&test_vcf_norm,$opts,in=>'norm.symbolic',fai=>'norm.symbolic',out=>'norm.symbolic.1.out',args=>'--old-rec-tag ORI');
+run_test(\&test_vcf_norm,$opts,in=>'norm.symbolic.2',fai=>'norm.symbolic',out=>'norm.symbolic.2.out',args=>'--old-rec-tag ORI');
 run_test(\&test_vcf_norm,$opts,in=>'norm.right-align',fai=>'norm.right-align',out=>'norm.right-align.1.out',args=>'--old-rec-tag ORI');
 run_test(\&test_vcf_norm,$opts,in=>'norm.right-align',fai=>'norm.right-align',out=>'norm.right-align.2.out',args=>'--old-rec-tag ORI -g {PATH}/norm.right-align.gff');
 run_test(\&test_vcf_view,$opts,in=>'weird-chr-names',out=>'weird-chr-names.1.out',args=>'',reg=>'-r 1');

vcfnorm.c

@@ -559,6 +559,7 @@ static int realign(args_t *args, bcf1_t *line)
     hts_expand0(kstring_t,line->n_allele,args->ntmp_del,args->tmp_del);
     kstring_t *als = args->tmp_als;
     kstring_t *del = args->tmp_del;
+    int symbolic_alts = 1;
     for (i=0; i<line->n_allele; i++)
     {
         del[i].l = 0;
@@ -576,12 +577,13 @@ static int realign(args_t *args, bcf1_t *line)
                 replace_iupac_codes(ref,nref);  // any non-ACGT character in fasta ref is replaced with N
                 als[0].l = 0;
                 kputs(ref, &als[0]);
-                als[i].l = 0;
-                kputsn(ref,1,&als[i]);
-                kputs(line->d.allele[i],&del[i]);
-                continue;
             }
+            als[i].l = 0;
+            kputsn(als[0].s,1,&als[i]);
+            kputs(line->d.allele[i],&del[i]);
+            continue;
         }
+        if ( i>0 ) symbolic_alts = 0;
         if ( line->d.allele[i][0]=='*' ) return ERR_SPANNING_DELETION;  // spanning deletion
         if ( has_non_acgtn(line->d.allele[i],line->shared.l) )
         {
@@ -610,8 +612,15 @@ static int realign(args_t *args, bcf1_t *line)
     else
         new_pos = realign_right(args, line);
 
-    // Have the alleles changed?
-    als[0].s[ als[0].l ] = 0;  // in order for strcmp to work
+    // Have the alleles changed? Consider <DEL> could have expanded the REF allele. In that
+    // case it must be trimmed, however the new REF length must reflect the entire length.
+    als[0].s[ als[0].l ] = 0;   // for strcmp to work
+    int new_reflen = strlen(als[0].s);
+    if ( symbolic_alts )
+    {
+        als[0].l = 1;
+        als[0].s[ als[0].l ] = 0;
+    }
     if ( new_pos==line->pos && !strcasecmp(line->d.allele[0],als[0].s) ) return ERR_OK;
 
     set_old_rec_tag(args, line, line, 0);
@@ -629,7 +638,6 @@ static int realign(args_t *args, bcf1_t *line)
     args->nchanged++;
 
     // Update INFO/END if necessary
-    int new_reflen = strlen(line->d.allele[0]);
     if ( (new_pos!=line->pos || reflen!=new_reflen) && bcf_get_info_int32(args->hdr, line, "END", &args->int32_arr, &args->nint32_arr)==1 )
     {
         // bcf_update_alleles_str() messed up rlen because line->pos changed. This will be fixed by bcf_update_info_int32()