When first base of the reference genome is deleted,

the VCF record has POS=1 and the first REF base cannot precede the event. Fixes #1330
samtools · Oct 27, 2020 · 9c15769 · 9c15769
1 parent feb1ffc
commit 9c15769
Show file tree

Hide file tree

Showing 5 changed files with 18 additions and 2 deletions.
diff --git a/consensus.c b/consensus.c
@@ -580,7 +580,14 @@ static void apply_variant(args_t *args, bcf1_t *rec)
     int trim_beg = 0;
     int var_type = bcf_get_variant_type(rec,ialt);
     int var_len  = rec->d.var[ialt].n;
-    if ( var_type & VCF_INDEL ) trim_beg = 1;
+    if ( var_type & VCF_INDEL )
+    {
+        // normally indel starts one base after, but not if the first base of the fa reference is deleted
+        if ( rec->d.allele[0][0] == rec->d.allele[ialt][0] )
+            trim_beg = 1;
+        else
+            trim_beg = 0;
+    }
     else if ( (var_type & VCF_OTHER) && !strcasecmp(rec->d.allele[ialt],"<DEL>") )
     {
         trim_beg = 1;
@@ -698,7 +705,6 @@ static void apply_variant(args_t *args, bcf1_t *rec)
     if ( len_diff <= 0 )
     {
         // deletion or same size event
-
         assert( args->fa_buf.l >= idx+rec->rlen );
         args->prev_base = args->fa_buf.s[idx+rec->rlen-1];
         args->prev_base_pos = rec->pos + rec->rlen - 1;

diff --git a/test/consensus.14.fa b/test/consensus.14.fa
@@ -0,0 +1,2 @@
+>1
+GACT
diff --git a/test/consensus.14.out b/test/consensus.14.out
@@ -0,0 +1,2 @@
+>1
+ACT
diff --git a/test/consensus.14.vcf b/test/consensus.14.vcf
@@ -0,0 +1,5 @@
+##fileformat=VCFv4.2
+##reference=genome.fa
+##contig=<ID=1,length=249250621>
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	.
+1	1	.	GA	A	.	.	.
diff --git a/test/test.pl b/test/test.pl
@@ -602,6 +602,7 @@
 test_vcf_consensus($opts,in=>'consensus.11',out=>'consensus.11.2.out',fa=>'consensus.11.fa',args=>q[-s smpl -a N]);
 test_vcf_consensus($opts,in=>'consensus.12',out=>'consensus.12.out',fa=>'consensus.12.fa',args=>'');
 test_vcf_consensus($opts,in=>'consensus.13',out=>'consensus.13.out',fa=>'consensus.13.fa',args=>'');
+test_vcf_consensus($opts,in=>'consensus.14',out=>'consensus.14.out',fa=>'consensus.14.fa',args=>'');
 test_mpileup($opts,in=>[qw(mpileup.1 mpileup.2 mpileup.3)],out=>'mpileup/mpileup.1.out',args=>q[-r17:100-150],test_list=>1);
 test_mpileup($opts,in=>[qw(mpileup.1 mpileup.2 mpileup.3)],out=>'mpileup/mpileup.2.out',args=>q[-a DP,DV -r17:100-600]); # test files from samtools mpileup test suite
 test_mpileup($opts,in=>[qw(mpileup.1)],out=>'mpileup/mpileup.3.out',args=>q[-B --ff 0x14 -r17:1050-1060]); # test file converted to vcf from samtools mpileup test suite