New -m snp-ins-del switch to merge SNVs, insertions and deletions separately

pd3 · jkbonfield · commit 6c66ae3fd79a · 2022-08-11T15:08:07.000+01:00
Resolves #1704 and depends on samtools/htslib#1467 to be merged
diff --git a/NEWS b/NEWS
@@ -24,6 +24,10 @@ Changes affecting specific commands:
       requested, which would result in not transferring the annotation from a tab-delimited
       file (#1733)
 
+* bcftools merge
+
+    - New `-m snp-ins-del` switch to merge SNVs, insertions and deletions separately (#1704)
+
 * bcftools mpileup
 
     - New NMBZ annotation for Mann-Whitney U-z test on number of mismatches within
diff --git a/doc/bcftools.txt b/doc/bcftools.txt
@@ -1902,15 +1902,16 @@ For "vertical" merge take a look at *<<concat,bcftools concat>>* or *<<norm,bcft
     maximum number of alternate alleles that can be included in the PL tag. The default value
     is 0 which disables the feature and outputs values for all alternate alleles.
 
-*-m, --merge* 'snps'|'indels'|'both'|'all'|'none'|'id'::
+*-m, --merge* 'snps'|'indels'|'both'|'snp-ins-del'|'all'|'none'|'id'::
     The option controls what types of multiallelic records can be created:
 ----
--m none   ..  no new multiallelics, output multiple records instead
--m snps   ..  allow multiallelic SNP records
--m indels ..  allow multiallelic indel records
--m both   ..  both SNP and indel records can be multiallelic
--m all    ..  SNP records can be merged with indel records
--m id     ..  merge by ID
+-m none        ..  no new multiallelics, output multiple records instead
+-m snps        ..  allow multiallelic SNP records
+-m indels      ..  allow multiallelic indel records
+-m both        ..  both SNP and indel records can be multiallelic
+-m all         ..  SNP records can be merged with indel records
+-m snp-ins-del ..  allow multiallelic SNVs, insertions, deletions, but don't mix them
+-m id          ..  merge by ID
 ----
 
 *--no-index*::
diff --git a/plugins/smpl-stats.c b/plugins/smpl-stats.c
@@ -1,19 +1,19 @@
 /* The MIT License
 
-   Copyright (c) 2018-2021 Genome Research Ltd.
+   Copyright (c) 2018-2022 Genome Research Ltd.
 
    Author: Petr Danecek <pd3@sanger.ac.uk>
-   
+
    Permission is hereby granted, free of charge, to any person obtaining a copy
    of this software and associated documentation files (the "Software"), to deal
    in the Software without restriction, including without limitation the rights
    to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
    copies of the Software, and to permit persons to whom the Software is
    furnished to do so, subject to the following conditions:
-   
+
    The above copyright notice and this permission notice shall be included in
    all copies or substantial portions of the Software.
-   
+
    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
    IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
    FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
@@ -92,7 +92,7 @@ const char *about(void)
 
 static const char *usage_text(void)
 {
-    return 
+    return
         "\n"
         "About: Calculates basic per-sample stats. Use curly brackets to scan a range of values simultaneously\n"
         "Usage: bcftools +smpl-stats [Plugin Options]\n"
@@ -147,7 +147,7 @@ static void parse_filters(args_t *args)
         }
         if ( !expanded ) break;
     }
-    
+
     fprintf(stderr,"Collecting data for %d filtering expressions\n", args->nflt_str);
 }
 
@@ -184,9 +184,9 @@ static void init_data(args_t *args)
             // replace tab's with spaces so that the output stays parsable
             char *tmp = args->filters[i].expr;
             while ( *tmp )
-            { 
-                if ( *tmp=='\t' ) *tmp = ' '; 
-                tmp++; 
+            {
+                if ( *tmp=='\t' ) *tmp = ' ';
+                tmp++;
             }
         }
     }
@@ -339,7 +339,7 @@ static void process_record(args_t *args, bcf1_t *rec, flt_stats_t *flt)
     int ngt = bcf_get_genotypes(args->hdr, rec, &args->gt_arr, &args->mgt_arr);
     if ( ngt<0 ) return;
     int ngt1 = ngt / rec->n_sample;
-    
+
 
     // For ts/tv: numeric code of the reference allele, -1 for insertions
     int ref = !rec->d.allele[0][1] ? bcf_acgt2int(*rec->d.allele[0]) : -1;
@@ -381,7 +381,7 @@ static void process_record(args_t *args, bcf1_t *rec, flt_stats_t *flt)
             has_nonref = 1;
         }
         if ( !has_nonref ) continue; // only ref or * in this genotype
-        
+
         stats->nnon_ref++;
 
         // Calculate ts/tv, count SNPs, indels. It does the right thing and handles also HetAA genotypes
@@ -395,8 +395,7 @@ static void process_record(args_t *args, bcf1_t *rec, flt_stats_t *flt)
 
                 if ( args->ac[als[j]]==1 ) { stats->nsingleton++; site_singleton = 1; }
 
-                int var_type = bcf_get_variant_type(rec, als[j]);
-                if ( var_type==VCF_SNP || var_type==VCF_MNP )
+                if ( bcf_has_variant_type(rec, als[j], VCF_SNP|VCF_MNP, subset) )
                 {
                     int k = 0;
                     while ( rec->d.allele[0][k] && rec->d.allele[als[j]][k] )
@@ -411,7 +410,7 @@ static void process_record(args_t *args, bcf1_t *rec, flt_stats_t *flt)
                         k++;
                     }
                 }
-                else if ( var_type==VCF_INDEL ) has_indel = 1;
+                else if ( bcf_has_variant_type(rec, als[j], VCF_INDEL, exact) ) has_indel = 1;
             }
             if ( has_ts ) { stats->nts++; site_has_ts = 1; }
             if ( has_tv ) { stats->ntv++; site_has_tv = 1; }
@@ -446,7 +445,7 @@ int run(int argc, char **argv)
     int c, i;
     while ((c = getopt_long(argc, argv, "o:s:i:e:r:R:t:T:",loptions,NULL)) >= 0)
     {
-        switch (c) 
+        switch (c)
         {
             case 'e':
                 if ( args->filter_str ) error("Error: only one -i or -e expression can be given, and they cannot be combined\n");
diff --git a/test/merge.10.1.out b/test/merge.10.1.out
@@ -0,0 +1,13 @@
+##fileformat=VCFv4.3
+##FILTER=<ID=PASS,Description="All filters passed">
+##contig=<ID=1,assembly=b37,length=249250621>
+##reference=file:ref.fa
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
+1	3000000	.	C	CCG	.	.	.
+1	3000000	.	C	CAA	.	.	.
+1	3000150	.	CC	C	.	.	.
+1	3000150	.	C	CTT	.	.	.
+1	3000152	.	C	A	.	.	.
+1	3000152	.	C	CA	.	.	.
+1	3000154	.	C	A	.	.	.
+1	3000154	.	C	T	.	.	.
diff --git a/test/merge.10.2.out b/test/merge.10.2.out
@@ -0,0 +1,10 @@
+##fileformat=VCFv4.3
+##FILTER=<ID=PASS,Description="All filters passed">
+##contig=<ID=1,assembly=b37,length=249250621>
+##reference=file:ref.fa
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
+1	3000000	.	C	CCG,CAA	.	.	.
+1	3000150	.	CC	C,CTTC	.	.	.
+1	3000152	.	C	A	.	.	.
+1	3000152	.	C	CA	.	.	.
+1	3000154	.	C	A,T	.	.	.
diff --git a/test/merge.10.3.out b/test/merge.10.3.out
@@ -0,0 +1,11 @@
+##fileformat=VCFv4.3
+##FILTER=<ID=PASS,Description="All filters passed">
+##contig=<ID=1,assembly=b37,length=249250621>
+##reference=file:ref.fa
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
+1	3000000	.	C	CCG,CAA	.	.	.
+1	3000150	.	CC	C,CTTC	.	.	.
+1	3000152	.	C	A	.	.	.
+1	3000152	.	C	CA	.	.	.
+1	3000154	.	C	A	.	.	.
+1	3000154	.	C	T	.	.	.
diff --git a/test/merge.10.a.vcf b/test/merge.10.a.vcf
@@ -0,0 +1,8 @@
+##fileformat=VCFv4.3
+##contig=<ID=1,assembly=b37,length=249250621>
+##reference=file:ref.fa
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
+1	3000000	.	C	CCG	.	.	.
+1	3000150	.	CC	C	.	.	.
+1	3000152	.	C	A	.	.	.
+1	3000154	.	C	A	.	.	.
diff --git a/test/merge.10.b.vcf b/test/merge.10.b.vcf
@@ -0,0 +1,8 @@
+##fileformat=VCFv4.3
+##contig=<ID=1,assembly=b37,length=249250621>
+##reference=file:ref.fa
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
+1	3000000	.	C	CAA	.	.	.
+1	3000150	.	C	CTT	.	.	.
+1	3000152	.	C	CA	.	.	.
+1	3000154	.	C	T	.	.	.
diff --git a/test/test.pl b/test/test.pl
@@ -97,6 +97,9 @@
 run_test(\&test_vcf_merge,$opts,in=>['merge.8.a','merge.8.b'],out=>'merge.8.out',args=>'-i AN:sum,AC:sum');
 run_test(\&test_vcf_merge,$opts,in=>['merge.9.a','merge.9.b'],out=>'merge.9.1.out',args=>'');
 run_test(\&test_vcf_merge,$opts,in=>['merge.9.a','merge.9.b'],out=>'merge.9.2.out',args=>'-i AN:sum,AC:sum');
+run_test(\&test_vcf_merge,$opts,in=>['merge.10.a','merge.10.b'],out=>'merge.10.1.out',args=>'-m none');
+run_test(\&test_vcf_merge,$opts,in=>['merge.10.a','merge.10.b'],out=>'merge.10.2.out',args=>'-m both');
+run_test(\&test_vcf_merge,$opts,in=>['merge.10.a','merge.10.b'],out=>'merge.10.3.out',args=>'-m snp-ins-del');
 # run_test(\&test_vcf_merge_big,$opts,in=>'merge_big.1',out=>'merge_big.1.1',nsmpl=>79000,nfiles=>79,nalts=>486,args=>'');   # commented out for speed
 run_test(\&test_vcf_query,$opts,in=>'query.string',out=>'query.string.1.out',args=>q[-f '%CHROM\\t%POS\\t%CLNREVSTAT\\n' -i'CLNREVSTAT="criteria_provided,_conflicting_interpretations"']);
 run_test(\&test_vcf_query,$opts,in=>'query.string',out=>'query.string.1.out',args=>q[-f '%CHROM\\t%POS\\t%CLNREVSTAT\\n' -i'CLNREVSTAT="criteria_provided" || CLNREVSTAT="_conflicting_interpretations"']);
diff --git a/vcfmerge.c b/vcfmerge.c
@@ -43,6 +43,8 @@ THE SOFTWARE.  */
 
 #define DBG 0
 
+#define COLLAPSE_SNP_INS_DEL (1<<10)
+
 #include <htslib/khash.h>
 KHASH_MAP_INIT_STR(strdict, int)
 typedef khash_t(strdict) strdict_t;
@@ -2517,7 +2519,16 @@ static inline int is_gvcf_block(bcf1_t *line)
     }
     return 0;
 }
-static const int snp_mask = (VCF_SNP<<2)|(VCF_MNP<<2), indel_mask = VCF_INDEL<<2, ref_mask = 2;
+
+// Lines can come with any combination of variant types. We use a subset of types defined in vcf.h
+// but shift by two bits to account for VCF_REF defined as 0 (design flaw in vcf.h, my fault) and
+// to accommodate for VCF_GVCF_REF defined below
+static const int
+    snp_mask = (VCF_SNP<<2)|(VCF_MNP<<2),
+    indel_mask = VCF_INDEL<<2,
+    ins_mask = VCF_INS<<2,
+    del_mask = VCF_DEL<<2,
+    ref_mask = 2;
 
 /*
     Check incoming lines for new gVCF blocks, set pointer to the current source
@@ -2743,6 +2754,11 @@ int can_merge(args_t *args)
             else
             {
                 int var_type = bcf_get_variant_types(line);
+                if ( args->collapse==COLLAPSE_SNP_INS_DEL )
+                {
+                    // need to distinguish between ins and del so strip the VCF_INDEL flag
+                    var_type &= ~VCF_INDEL;
+                }
                 maux->var_types |= var_type ? var_type<<2 : 2;
 
                 // for the `-m none -g` mode
@@ -2812,7 +2828,7 @@ int can_merge(args_t *args)
                     //  - SNPs+SNPs+MNPs+REF if -m both,snps
                     //  - indels+indels+REF  if -m both,indels, REF only if SNPs are not present
                     //  - SNPs come first
-                    if ( line_type & indel_mask )
+                    if ( line_type & (indel_mask|ins_mask|del_mask) )
                     {
                         if ( !(line_type&snp_mask) && maux->var_types&snp_mask ) continue;  // SNPs come first
                         if ( args->do_gvcf && maux->var_types&ref_mask ) continue;  // never merge indels with gVCF blocks
@@ -2898,16 +2914,22 @@ void stage_line(args_t *args)
                 int line_type = bcf_get_variant_types(buf->lines[j]);
                 if ( maux->var_types&snp_mask && line_type&VCF_SNP && (args->collapse&COLLAPSE_SNPS) ) break;
                 if ( maux->var_types&indel_mask && line_type&VCF_INDEL && (args->collapse&COLLAPSE_INDELS) ) break;
+                if ( maux->var_types&ins_mask && line_type&VCF_INS && (args->collapse&COLLAPSE_SNP_INS_DEL) ) break;
+                if ( maux->var_types&del_mask && line_type&VCF_DEL && (args->collapse&COLLAPSE_SNP_INS_DEL) ) break;
                 if ( line_type==VCF_REF )
                 {
                     if ( maux->var_types&snp_mask && (args->collapse&COLLAPSE_SNPS) ) break;
                     if ( maux->var_types&indel_mask && (args->collapse&COLLAPSE_INDELS) ) break;
+                    if ( maux->var_types&ins_mask && (args->collapse&COLLAPSE_SNP_INS_DEL) ) break;
+                    if ( maux->var_types&del_mask && (args->collapse&COLLAPSE_SNP_INS_DEL) ) break;
                     if ( maux->var_types&ref_mask ) break;
                 }
                 else if ( maux->var_types&ref_mask )
                 {
                     if ( line_type&snp_mask && (args->collapse&COLLAPSE_SNPS) ) break;
                     if ( line_type&indel_mask && (args->collapse&COLLAPSE_INDELS) ) break;
+                    if ( line_type&ins_mask && (args->collapse&COLLAPSE_SNP_INS_DEL) ) break;
+                    if ( line_type&del_mask && (args->collapse&COLLAPSE_SNP_INS_DEL) ) break;
                 }
             }
         }
@@ -3125,7 +3147,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           EXPERIMENTAL: 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|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|none|id>, see man page for details [both]\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");
     fprintf(stderr, "    -o, --output FILE                 Write output to a file [standard output]\n");
@@ -3229,6 +3251,7 @@ int main_vcfmerge(int argc, char *argv[])
                 else if ( !strcmp(optarg,"any") ) args->collapse |= COLLAPSE_ANY;
                 else if ( !strcmp(optarg,"all") ) args->collapse |= COLLAPSE_ANY;
                 else if ( !strcmp(optarg,"none") ) args->collapse = COLLAPSE_NONE;
+                else if ( !strcmp(optarg,"snp-ins-del") ) args->collapse = COLLAPSE_SNP_INS_DEL;
                 else if ( !strcmp(optarg,"id") ) { args->collapse = COLLAPSE_NONE; args->merge_by_id = 1; }
                 else error("The -m type \"%s\" is not recognised.\n", optarg);
                 break;
diff --git a/vcfstats.c b/vcfstats.c
@@ -1,6 +1,6 @@
 /*  vcfstats.c -- Produces stats which can be plotted using plot-vcfstats.
 
-    Copyright (C) 2012-2021 Genome Research Ltd.
+    Copyright (C) 2012-2022 Genome Research Ltd.
 
     Author: Petr Danecek <pd3@sanger.ac.uk>
 
@@ -706,7 +706,7 @@ static void do_indel_stats(args_t *args, stats_t *stats, bcf_sr_t *reader)
     for (i=1; i<line->n_allele; i++)
     {
         if ( args->first_allele_only && i>1 ) break;
-        if ( bcf_get_variant_type(line,i)!=VCF_INDEL ) continue;
+        if ( !bcf_has_variant_type(line,i,VCF_INDEL,exact) ) continue;
         int len = line->d.var[i].n;
 
         #if IRC_STATS
