Careful when merging symbolic alleles

Don't merge symbolic alleles unless they have not just the same type, eg. <DEL>,
but also length, i.e the INFO/END coordinate

Fixes #2362

Author: pd3

NEWS

@@ -42,6 +42,9 @@ Changes affecting specific commands:
     - Make `--merge both` work with indel-only records; for example, the multiallelic
       site G>GT,T should be merged with G>GT (#2339)
 
+    - 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)
+
 * bcftools norm
 
     - Print the number of removed duplicate sites in the final statistics (#2346)

test/merge.symbolic.1.1.out

@@ -0,0 +1,8 @@
+##fileformat=VCFv4.1
+##FILTER=<ID=PASS,Description="All filters passed">
+##contig=<ID=chr1,length=248956422>
+##INFO=<ID=END,Number=1,Type=Integer,Description="End position of structural variation">
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	Sample	Other
+chr1	147022730	SV1	N	<DEL>	.	.	END=147593064	GT	0/1	./.
+chr1	147022730	SV2	N	<DEL>	.	.	END=148013144	GT	./.	1/1

test/merge.symbolic.1.a.vcf

@@ -0,0 +1,6 @@
+##fileformat=VCFv4.1
+##contig=<ID=chr1,length=248956422>
+##INFO=<ID=END,Number=1,Type=Integer,Description="End position of structural variation">
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	Sample
+chr1	147022730	SV1	N	<DEL>	.	.	END=147593064	GT	0/1

test/merge.symbolic.1.b.vcf

@@ -0,0 +1,6 @@
+##fileformat=VCFv4.1
+##contig=<ID=chr1,length=248956422>
+##INFO=<ID=END,Number=1,Type=Integer,Description="End position of structural variation">
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	Other
+chr1	147022730	SV2	N	<DEL>	.	.	END=148013144	GT	1/1

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.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');
 run_test(\&test_vcf_merge,$opts,in=>['merge.12.a','merge.12.b'],out=>'merge.12.1.out',args=>'--merge none');

vcfmerge.c

@@ -132,6 +132,7 @@ typedef struct
     int mrec;       // allocated size of buf
     maux1_t *rec;   // buffer to keep reader's lines
     bcf1_t **lines; // source buffer: either gvcf or readers' buffer
+    bcf_hdr_t *hdr; // this reader's header
     int var_types;  // reader's variant types in the active [beg,end] window
 }
 buffer_t;
@@ -871,7 +872,10 @@ maux_t *maux_init(args_t *args)
     ma->smpl_nGsize = (int*) malloc(n_smpl*sizeof(int));
     ma->buf = (buffer_t*) calloc(ma->n,sizeof(buffer_t));
     for (i=0; i<ma->n; i++)
+    {
         ma->buf[i].rid = -1;
+        ma->buf[i].hdr = files->readers[i].header;
+    }
     ma->str = (kstring_t*) calloc(n_smpl,sizeof(kstring_t));
     if ( args->local_alleles )
     {
@@ -2925,23 +2929,48 @@ static const int
     indel_mask = (VCF_INDEL<<1),
     ins_mask   = VCF_INS<<1,
     del_mask   = VCF_DEL<<1,
-    ref_mask   = 1;
+    ref_mask   = 1,
+    other_mask = VCF_OTHER<<1;
+
+typedef struct
+{
+    int types,      // selected types, see the *_mask(s) above
+        end;        // if symbolic allele is involved, the END coordinate of the first record
+    bcf1_t *rec;    // the first record selected
+}
+selected_t;
 
 // Can these types be merged given the -m settings? Despite the function's name, its focus is on
 // excluding incompatible records, there will be a finer matching later in stage_line()
-static inline int types_compatible(args_t *args, int selected_types, buffer_t *buf, int irec)
+static inline int types_compatible(args_t *args, selected_t *selected, buffer_t *buf, int irec)
 {
     int k;
     maux_t *maux = args->maux;
     bcf1_t *rec = buf->lines[irec];
     int rec_types = buf->rec[irec].var_types;
 
-    assert( selected_types );   // this is trivially true, set in can_merge()
+    int end = -1;
+    if ( rec_types&other_mask )
+    {
+        int32_t *itmp = NULL, nitmp = 0;
+        bcf_get_info_int32(buf->hdr,rec,"END",&itmp,&nitmp);
+        end = nitmp==1 ? itmp[0] : -1;
+        free(itmp);
+    }
+
+    // First time here?
+    if ( !selected->types )
+    {
+        selected->end = end;
+        selected->rec = rec;
+        selected->types = rec_types;
+        return 1;
+    }
 
     if ( args->collapse & COLLAPSE_ANY ) return 1;  // can merge anything with anything
 
     // REF and gVCF_REF with no other alleles present can be merged with anything
-    if ( (selected_types&ref_mask) && !(selected_types&(~ref_mask)) ) return 1;
+    if ( (selected->types&ref_mask) && !(selected->types&(~ref_mask)) ) return 1;
     if ( (rec_types&ref_mask) && !(rec_types&(~ref_mask)) ) return 1;
 
     if ( args->collapse!=COLLAPSE_NONE )
@@ -2952,23 +2981,23 @@ static inline int types_compatible(args_t *args, int selected_types, buffer_t *b
         //  - rec has indel, we already have an indel, and -m both,indels,snp-ins-del
         if ( args->collapse&(COLLAPSE_SNPS|COLLAPSE_SNP_INS_DEL) )
         {
-            if ( (rec_types&snp_mask) && (selected_types&snp_mask) ) return 1;
+            if ( (rec_types&snp_mask) && (selected->types&snp_mask) ) return 1;
         }
         if ( args->collapse&COLLAPSE_INDELS )
         {
-            if ( (rec_types&indel_mask) && (selected_types&indel_mask) ) return 1;
+            if ( (rec_types&indel_mask) && (selected->types&indel_mask) ) return 1;
         }
         if ( args->collapse&COLLAPSE_SNP_INS_DEL )
         {
-            if ( (rec_types&ins_mask) && (selected_types&ins_mask) ) return 1;
-            if ( (rec_types&del_mask) && (selected_types&del_mask) ) return 1;
+            if ( (rec_types&ins_mask) && (selected->types&ins_mask) ) return 1;
+            if ( (rec_types&del_mask) && (selected->types&del_mask) ) return 1;
         }
         // Whatever is left, allow to match if the alleles match exactly
     }
 
     // The -m none mode or exact matching requested
     // Simple test first: are the variants of the same type?
-    int x = selected_types;
+    int x = selected->types;
     int y = rec_types;
     if ( !(x&y) ) return 0;                 // no matching type
     if ( (x&y)!=x && (x&y)!=y ) return 0;   // not a subset
@@ -2980,6 +3009,13 @@ static inline int types_compatible(args_t *args, int selected_types, buffer_t *b
         if ( vcmp_find_allele(args->vcmp,maux->als+1,maux->nals-1,rec->d.allele[k])>=0 ) break;
     }
     if ( k==rec->n_allele ) return 0;   // this record has a new allele rec->d.allele[k]
+
+    if ( selected->types&other_mask && rec_types&other_mask )
+    {
+        // both records have symbolic alleles and the alleles are the same
+        if ( selected->end!=end ) return 0;
+    }
+
     return 1;   // all alleles in rec are also in the records selected thus far, perhaps save for gVCF_REF
 }
 
@@ -3106,7 +3142,7 @@ int can_merge(args_t *args)
     }
     if ( !ntodo ) return 0;
 
-    int selected_types = 0;
+    selected_t selected = {0,0,NULL};
 
     // In this loop we select from each reader compatible candidate lines.
     // (i.e. SNPs or indels). Go through all files and all lines at this
@@ -3121,7 +3157,7 @@ int can_merge(args_t *args)
             gaux[i].line->d.allele[0][0] = ref;
             gaux[i].line->pos = maux->pos;
             maux_update_alleles(args, i, buf->beg);
-            selected_types |= ref_mask;
+            selected.types |= ref_mask;
             continue;
         }
         for (j=buf->beg; j<buf->end; j++)
@@ -3136,7 +3172,7 @@ int can_merge(args_t *args)
             {
                 if ( strcmp(id,line->d.id) ) continue;      // matching by ID and it does not match the selected record
             }
-            else if ( selected_types && !types_compatible(args,selected_types,buf,j) ) continue;
+            else if ( !types_compatible(args,&selected,buf,j) ) continue;
 
             // This is not a good code. It makes the incorrect assumption of always having a SNP record available.
             // However, that is not always the case and prevents the merging of G>GT,T with G>GT (see test/merge.multiallelics.1.*.vcf).
@@ -3154,7 +3190,7 @@ int can_merge(args_t *args)
             //         ) continue;
             //  }
 
-            selected_types |= line_types;
+            selected.types |= line_types;
 
             buf->rec[j].skip = 0;   // the j-th record from i-th reader can be included. Final decision will be made in stage_line
             maux_update_alleles(args, i, j);