Add new -m exact merging mode. Resolves #2498

Author: pd3

NEWS

@@ -25,6 +25,9 @@ Changes affecting specific commands:
     - Fixes in memory access with localized alleles, it could result in printing
       printing trailing garbage in localized FORMAT fields (#2520)
 
+    - Add new `-m exact` mode to merge only records with the exact same alleles
+      (#2498)
+
 * bcftools mpileup
 
     - Remove unused experimental INFO/MIN_PL_SUM annotation

doc/bcftools.txt

@@ -2141,6 +2141,7 @@ For "vertical" merge take a look at *<<concat,bcftools concat>>* or *<<norm,bcft
     '\*' is appended, the unobserved allele '<*>' or '<NON_REF>' will be removed at variant sites;
     if two asterisks '**' are appended, the unobserved allele will be removed all sites.
 ----
+-m exact       ..  require the exact same alleles
 -m none        ..  no new multiallelics, output multiple records instead
 -m snps        ..  allow multiallelic SNP records
 -m indels      ..  allow multiallelic indel records

test/merge.exact.1.out

@@ -0,0 +1,7 @@
+##fileformat=VCFv4.2
+##FILTER=<ID=PASS,Description="All filters passed">
+##contig=<ID=1>
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
+1	100	.	C	T,G	.	.	.
+1	200	.	C	T,A,G	.	.	.
+1	300	.	C	T,T,G	.	.	.

test/merge.exact.2.out

@@ -0,0 +1,10 @@
+##fileformat=VCFv4.2
+##FILTER=<ID=PASS,Description="All filters passed">
+##contig=<ID=1>
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
+1	100	.	C	T	.	.	.
+1	100	.	C	T,G	.	.	.
+1	200	.	C	T,A	.	.	.
+1	200	.	C	T,G	.	.	.
+1	300	.	C	T,T	.	.	.
+1	300	.	C	T,G	.	.	.

test/merge.exact.a.vcf

@@ -0,0 +1,6 @@
+##fileformat=VCFv4.2
+##contig=<ID=1>
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
+1	100	.	C	T	.	.	.
+1	200	.	C	T,A	.	.	.
+1	300	.	C	T,T	.	.	.

test/merge.exact.b.vcf

@@ -0,0 +1,6 @@
+##fileformat=VCFv4.2
+##contig=<ID=1>
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
+1	100	.	C	T,G	.	.	.
+1	200	.	C	T,G	.	.	.
+1	300	.	C	T,G	.	.	.

test/test.pl

@@ -64,6 +64,8 @@
 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-phase.a','merge.broken-phase.b'],out=>'merge.broken-phase.1.out',args=>'');
+run_test(\&test_vcf_merge,$opts,in=>['merge.exact.a','merge.exact.b'],out=>'merge.exact.1.out',args=>'');
+run_test(\&test_vcf_merge,$opts,in=>['merge.exact.a','merge.exact.b'],out=>'merge.exact.2.out',args=>'-m exact');
 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');

vcfmerge.c

@@ -1,6 +1,6 @@
 /*  vcfmerge.c -- Merge multiple VCF/BCF files to create one multi-sample file.
 
-    Copyright (C) 2012-2025 Genome Research Ltd.
+    Copyright (C) 2012-2026 Genome Research Ltd.
 
     Author: Petr Danecek <pd3@sanger.ac.uk>
 
@@ -2980,13 +2980,18 @@ static inline int types_compatible(args_t *args, selected_t *selected, buffer_t
         return 1;
     }
 
+    if ( args->collapse==BCF_SR_PAIR_EXACT )
+    {
+        if ( maux->nals != rec->n_allele ) return 0;
+    }
+
     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 ( (rec_types&ref_mask) && !(rec_types&(~ref_mask)) ) return 1;
 
-    if ( args->collapse!=COLLAPSE_NONE )
+    if ( args->collapse!=COLLAPSE_NONE && args->collapse!=BCF_SR_PAIR_EXACT )
     {
         // If we are here, one the following modes must have been set: both,snps,indels,snp-ins-del
         // Include the new record if
@@ -3013,19 +3018,31 @@ static inline int types_compatible(args_t *args, selected_t *selected, buffer_t
     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
+    if ( (x&y)!=x && (x&y)!=y ) return 0;   // neither is a type subset of the other
+    if ( args->collapse==BCF_SR_PAIR_EXACT && ((x&y)!=x || (x&y)!=y) ) return 0;
 
     if ( vcmp_set_ref(args->vcmp,maux->als[0],rec->d.allele[0]) < 0 ) return 0;   // refs are not compatible
-    for (k=1; k<rec->n_allele; k++)
+    if ( args->collapse!=BCF_SR_PAIR_EXACT )
+    {
+        for (k=1; k<rec->n_allele; k++)
+        {
+            if ( bcf_has_variant_type(rec,k,VCF_REF) ) continue;    // this must be gVCF_REF (<*> or <NON_REF>)
+            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;   // none of the alleles in rec is present in maux
+    }
+    else
     {
-        if ( bcf_has_variant_type(rec,k,VCF_REF) ) continue;    // this must be gVCF_REF (<*> or <NON_REF>)
-        if ( vcmp_find_allele(args->vcmp,maux->als+1,maux->nals-1,rec->d.allele[k])>=0 ) break;
+        for (k=1; k<rec->n_allele; k++)
+        {
+            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;   // at least one of the alleles in rec is not present in maux
     }
-    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
+        // both records have symbolic alleles and the alleles are the same - does END match as well?
         if ( selected->end!=end ) return 0;
     }
 
@@ -3267,7 +3284,7 @@ void stage_line(args_t *args)
         if ( j>=buf->end )
         {
             // no matching allele found in this file
-            if ( args->collapse==COLLAPSE_NONE ) continue;  // exact matching requested, skip
+            if ( args->collapse==COLLAPSE_NONE || args->collapse==BCF_SR_PAIR_EXACT ) continue;  // exact matching requested, skip
 
             // choose something compatible to create a multiallelic site given the -m criteria
             for (j=buf->beg; j<buf->end; j++)
@@ -3546,7 +3563,7 @@ static void usage(void)
     fprintf(stderr, "    -i, --info-rules TAG:METHOD,..    Rules for merging INFO fields (method is one of sum,avg,min,max,join) or \"-\" to turn off the default [DP:sum,DP4:sum]\n");
     fprintf(stderr, "    -l, --file-list FILE              Read file names from the file\n");
     fprintf(stderr, "    -L, --local-alleles INT           If more than INT alt alleles are encountered, drop FMT/PL and output LAA+LPL instead; 0=unlimited [0]\n");
-    fprintf(stderr, "    -m, --merge STRING[*|**]          Allow multiallelic records for snps,indels,both,snp-ins-del,all,none,id,*,**; see man page for details [both]\n");
+    fprintf(stderr, "    -m, --merge STRING[*|**]          Allow multiallelic records for snps,indels,both,snp-ins-del,all,exact,none,id,*,**; see man page for details [both]\n");
     fprintf(stderr, "    -M, --missing-rules TAG:METHOD    Rules for replacing missing values in numeric vectors (.,0,max) when unknown allele <*> is not present [.]\n");
     fprintf(stderr, "        --no-index                    Merge unindexed files, the same chromosomal order is required and -r/-R are not allowed\n");
     fprintf(stderr, "        --no-version                  Do not append version and command line to the header\n");
@@ -3669,6 +3686,7 @@ int main_vcfmerge(int argc, char *argv[])
                 else if ( !strncmp(optarg,"all",len) ) args->collapse |= COLLAPSE_ANY;
                 else if ( !strncmp(optarg,"both",len) ) args->collapse |= COLLAPSE_BOTH;
                 else if ( !strncmp(optarg,"none",len) ) args->collapse = COLLAPSE_NONE;
+                else if ( !strncmp(optarg,"exact",len) ) args->collapse = BCF_SR_PAIR_EXACT;
                 else error("The -m type \"%s\" is not recognised.\n", optarg);
                 break;
             }