Add flexibility to FILTER column transfers

pd3 · pd3 · commit 43d954e96907 · 2023-11-29T12:19:03.000+01:00
Allow transfers within the same file, across files, and in combination: # transfer FILTER column to INFO tag NewTag; notice that the -a option is not present, therefore # B.bcf/FILTER is the source annotation bcftools annotate -c INFO/NewTag:=FILTER B.bcf # transfer FILTER column from A.bcf to INFO/NewTag in B.bcf; notice that the -a option is present, # therefore A.bcf/FILTER is the source annotation bcftools annotate -c INFO/NewTag:=FILTER -a A.bcf B.bcf # transfer B.bcf/FILTER column to INFO tag NewTag; notice that although the -a option is present, # the notation "./FILTER" tells the program to transfer the annotation locally and make B.bcf/FILTER # the source annotation bcftools annotate -c INFO/NewTag:=./FILTER -a A.bcf B.bcf # transfer B.bcf/FILTER column to INFO/NewTag, then A.bcf/FILTER to B.bcf/FILTER bcftools annotate -c INFO/NewTag:=./FILTER,FILTER -a A.bcf B.bcf # transfer A.bcf/FILTER to B.bcf/FILTER, then the new B.bcf/FILTER value to INFO/NewTag; notice # that due to the order of operation, FILTER and INFO/NewTag will be identical bcftools annotate -c FILTER,INFO/NewTag:=./FILTER -a A.bcf B.bcf Resolves #1187 (comment) #1187 (comment)
diff --git a/NEWS b/NEWS
@@ -21,6 +21,10 @@ Changes affecting specific commands:
     - Fix `bcftools annotate --mark-sites`, VCF sites overlapping regions in a BED file
       were not annotated (#1989)
 
+    - Add flexibility to FILTER column transfers and allow transfers within the same file,
+      across files, and in combination. For examples see
+        http://samtools.github.io/bcftools/howtos/annotate.html#transfer_filter_to_info
+
 * bcftools call
 
     - Output MIN_DP rather than MinDP in gVCF mode
diff --git a/doc/bcftools.txt b/doc/bcftools.txt
@@ -542,7 +542,7 @@ Add or remove annotations.
     # that INFO/END is already present in the VCF header.
     bcftools annotate -a annots.tab.gz  -c CHROM,POS,~ID,REF,ALT,INFO/END input.vcf
 
-    # For more examples see http://samtools.github.io/bcftools/howtos/annotate.html
+    # For (many) more examples see http://samtools.github.io/bcftools/howtos/annotate.html
 ----
 
 
diff --git a/test/annotate19.1.out b/test/annotate19.1.out
@@ -2,7 +2,8 @@
 ##FILTER=<ID=PASS,Description="All filters passed">
 ##contig=<ID=1,assembly=b37,length=249250621>
 ##reference=file:///lustre/scratch105/projects/g1k/ref/main_project/human_g1k_v37.fasta
+##FILTER=<ID=ori_filter,Description="Test">
 ##INFO=<ID=ID,Number=1,Type=String,Description="Transferred ID column">
 ##INFO=<ID=INFO_ID,Number=1,Type=String,Description="Test">
 #CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
-1	3000001	id_id;info_id	C	T	.	.	ID=id_id;INFO_ID=info_id
+1	3000001	id_id;info_id	C	T	.	ori_filter	ID=id_id;INFO_ID=info_id
diff --git a/test/annotate19.2.out b/test/annotate19.2.out
@@ -2,6 +2,7 @@
 ##FILTER=<ID=PASS,Description="All filters passed">
 ##contig=<ID=1,assembly=b37,length=249250621>
 ##reference=file:///lustre/scratch105/projects/g1k/ref/main_project/human_g1k_v37.fasta
+##FILTER=<ID=ori_filter,Description="Test">
 ##FILTER=<ID=filter_filter,Description="Test">
 ##INFO=<ID=FILTER,Number=1,Type=String,Description="Transferred FILTER column">
 ##INFO=<ID=INFO_FILTER,Number=1,Type=String,Description="Test">
diff --git a/test/annotate19.3.out b/test/annotate19.3.out
@@ -0,0 +1,9 @@
+##fileformat=VCFv4.1
+##FILTER=<ID=PASS,Description="All filters passed">
+##contig=<ID=1,assembly=b37,length=249250621>
+##reference=file:///lustre/scratch105/projects/g1k/ref/main_project/human_g1k_v37.fasta
+##FILTER=<ID=ori_filter,Description="Test">
+##INFO=<ID=FILTER,Number=1,Type=String,Description="Transferred FILTER column">
+##INFO=<ID=INFO_FILTER,Number=1,Type=String,Description="Test">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
+1	3000001	ori_id	C	T	.	ori_filter	FILTER=filter_filter;INFO_FILTER=info_filter
diff --git a/test/annotate19.4.out b/test/annotate19.4.out
@@ -0,0 +1,8 @@
+##fileformat=VCFv4.1
+##FILTER=<ID=PASS,Description="All filters passed">
+##contig=<ID=1,assembly=b37,length=249250621>
+##reference=file:///lustre/scratch105/projects/g1k/ref/main_project/human_g1k_v37.fasta
+##FILTER=<ID=ori_filter,Description="Test">
+##INFO=<ID=FILTER,Number=1,Type=String,Description="Transferred FILTER column">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
+1	3000001	ori_id	C	T	.	ori_filter	FILTER=ori_filter
diff --git a/test/annotate19.5.out b/test/annotate19.5.out
@@ -0,0 +1,8 @@
+##fileformat=VCFv4.1
+##FILTER=<ID=PASS,Description="All filters passed">
+##contig=<ID=1,assembly=b37,length=249250621>
+##reference=file:///lustre/scratch105/projects/g1k/ref/main_project/human_g1k_v37.fasta
+##FILTER=<ID=ori_filter,Description="Test">
+##INFO=<ID=FILTER,Number=1,Type=String,Description="Transferred FILTER column">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
+1	3000001	ori_id	C	T	.	ori_filter	FILTER=filter_filter
diff --git a/test/annotate19.6.out b/test/annotate19.6.out
@@ -0,0 +1,9 @@
+##fileformat=VCFv4.1
+##FILTER=<ID=PASS,Description="All filters passed">
+##contig=<ID=1,assembly=b37,length=249250621>
+##reference=file:///lustre/scratch105/projects/g1k/ref/main_project/human_g1k_v37.fasta
+##FILTER=<ID=ori_filter,Description="Test">
+##INFO=<ID=FILTER,Number=1,Type=String,Description="Transferred FILTER column">
+##FILTER=<ID=filter_filter,Description="Test">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
+1	3000001	ori_id	C	T	.	filter_filter	FILTER=ori_filter
diff --git a/test/annotate19.7.out b/test/annotate19.7.out
@@ -0,0 +1,9 @@
+##fileformat=VCFv4.1
+##FILTER=<ID=PASS,Description="All filters passed">
+##contig=<ID=1,assembly=b37,length=249250621>
+##reference=file:///lustre/scratch105/projects/g1k/ref/main_project/human_g1k_v37.fasta
+##FILTER=<ID=ori_filter,Description="Test">
+##FILTER=<ID=filter_filter,Description="Test">
+##INFO=<ID=FILTER,Number=1,Type=String,Description="Transferred FILTER column">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
+1	3000001	ori_id	C	T	.	filter_filter	FILTER=filter_filter
diff --git a/test/annotate19.dst.vcf b/test/annotate19.dst.vcf
@@ -1,5 +1,6 @@
 ##fileformat=VCFv4.1
 ##contig=<ID=1,assembly=b37,length=249250621>
 ##reference=file:///lustre/scratch105/projects/g1k/ref/main_project/human_g1k_v37.fasta
+##FILTER=<ID=ori_filter,Description="Test">
 #CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
-1	3000001	ori_id	C	T	.	.	.
+1	3000001	ori_id	C	T	.	ori_filter	.
diff --git a/test/test.pl b/test/test.pl
@@ -518,6 +518,12 @@
 run_test(\&test_vcf_annotate,$opts,in=>'annotate18.2',tab=>'annotate18.2',out=>'annotate18.2.out',args=>'-c CHROM,BEG,END,A,B,C,D,E -l A:sum,B:avg,C:min,D:max,E:append');
 run_test(\&test_vcf_annotate,$opts,in=>'annotate19.dst',vcf=>'annotate19.src',out=>'annotate19.1.out',args=>'-c INFO/ID:=ID,INFO/INFO_ID:=INFO/ID,ID,=ID:=INFO/ID');
 run_test(\&test_vcf_annotate,$opts,in=>'annotate19.dst',vcf=>'annotate19.src',out=>'annotate19.2.out',args=>'-c FILTER,INFO/FILTER:=FILTER,INFO/INFO_FILTER:=INFO/FILTER');
+run_test(\&test_vcf_annotate,$opts,in=>'annotate19.dst',vcf=>'annotate19.src',out=>'annotate19.3.out',args=>'-c INFO/FILTER:=FILTER,INFO/INFO_FILTER:=INFO/FILTER');
+run_test(\&test_vcf_annotate,$opts,in=>'annotate19.dst',                      out=>'annotate19.4.out',args=>'-c INFO/FILTER:=FILTER');
+run_test(\&test_vcf_annotate,$opts,in=>'annotate19.dst',vcf=>'annotate19.src',out=>'annotate19.5.out',args=>'-c INFO/FILTER:=FILTER');
+run_test(\&test_vcf_annotate,$opts,in=>'annotate19.dst',vcf=>'annotate19.src',out=>'annotate19.4.out',args=>'-c INFO/FILTER:=./FILTER');
+run_test(\&test_vcf_annotate,$opts,in=>'annotate19.dst',vcf=>'annotate19.src',out=>'annotate19.6.out',args=>'-c INFO/FILTER:=./FILTER,FILTER');
+run_test(\&test_vcf_annotate,$opts,in=>'annotate19.dst',vcf=>'annotate19.src',out=>'annotate19.7.out',args=>'-c FILTER,INFO/FILTER:=./FILTER');
 run_test(\&test_vcf_annotate,$opts,in=>'annotate20.dst',vcf=>'annotate20.src',out=>'annotate20.1.out',args=>'-c  FMT/GT');
 run_test(\&test_vcf_annotate,$opts,in=>'annotate20.dst',vcf=>'annotate20.src',out=>'annotate20.2.out',args=>'-c +FMT/GT');
 run_test(\&test_vcf_annotate,$opts,in=>'annotate20.dst',vcf=>'annotate20.src',out=>'annotate20.3.out',args=>'-c -FMT/GT');
diff --git a/vcfannotate.c b/vcfannotate.c
@@ -170,6 +170,7 @@ typedef struct _args_t
     int argc, drop_header, record_cmd_line, tgts_is_vcf, mark_sites_logic, force, single_overlaps;
     int columns_is_file, has_append_mode, pair_logic;
     dbuf_t *header_lines;
+    bcf1_t *current_rec;    // current record for local setters
 }
 args_t;
 
@@ -517,8 +518,35 @@ static int vcf_getter_id2str(args_t *args, bcf1_t *rec, annot_col_t *col, void *
     *mptr = len+1;
     return len;
 }
+static int vcf_getter_filter2str_local(args_t *args, bcf1_t *rec, annot_col_t *col, void **ptr, int *mptr)
+{
+    rec = args->current_rec;
+    if ( !(rec->unpacked & BCF_UN_FLT) ) bcf_unpack(rec, BCF_UN_FLT);
+
+    kstring_t str;
+    str.s = *((char**)ptr);
+    str.m = *mptr;
+    str.l = 0;
+
+    int i;
+    if ( rec->d.n_flt )
+    {
+        for (i=0; i<rec->d.n_flt; i++)
+        {
+            if (i) kputc(';', &str);
+            kputs(bcf_hdr_int2id(args->hdr_out,BCF_DT_ID,rec->d.flt[i]), &str);
+        }
+    }
+    else kputc('.', &str);
+
+    *((char**)ptr) = str.s;
+    *mptr = str.m;
+    return str.l;
+}
 static int vcf_getter_filter2str(args_t *args, bcf1_t *rec, annot_col_t *col, void **ptr, int *mptr)
 {
+    if ( !(rec->unpacked & BCF_UN_FLT) ) bcf_unpack(rec, BCF_UN_FLT);
+
     kstring_t str;
     str.s = *((char**)ptr);
     str.m = *mptr;
@@ -2290,10 +2318,30 @@ static void init_columns(args_t *args)
             if ( bcf_hdr_id2type(args->tgts_hdr,BCF_HL_INFO,hdr_id)!=BCF_HT_STR )
                 error("Only Type=String tags can be used to annotate the ID column\n");
         }
-        else if ( (ptr=strstr(str.s,":=")) && !args->targets_fname )
+        else if ( (ptr=strstr(str.s,":=")) && (!args->targets_fname || !strncasecmp(ptr+2,"./",2)) )
         {
             *ptr = 0;
-            rename_annots_push(args,ptr+2,str.s);
+            if ( !strncasecmp(str.s,"INFO/",5) && (!strcasecmp(ptr+2,"FILTER") || !strcasecmp(ptr+2,"./FILTER")) )
+            {
+                // -a not present and transferring filter, needs to be a local transfer
+                args->ncols++; args->cols = (annot_col_t*) realloc(args->cols,sizeof(annot_col_t)*args->ncols);
+                annot_col_t *col = &args->cols[args->ncols-1];
+                memset(col,0,sizeof(*col));
+                col->icol = icol;
+                col->replace = replace;
+                col->setter = vcf_setter_info_str;
+                col->getter = vcf_getter_filter2str_local;
+                col->hdr_key_src = strdup(ptr+2);
+                col->hdr_key_dst = strdup(str.s+5);
+                tmp.l = 0;
+                ksprintf(&tmp,"##INFO=<ID=%s,Number=1,Type=String,Description=\"Transferred FILTER column\">",col->hdr_key_dst);
+                bcf_hdr_append(args->hdr_out, tmp.s);
+                if (bcf_hdr_sync(args->hdr_out) < 0) error_errno("[%s] Failed to update header", __func__);
+                int hdr_id = bcf_hdr_id2int(args->hdr_out, BCF_DT_ID, col->hdr_key_dst);
+                col->number = bcf_hdr_id2length(args->hdr_out,BCF_HL_INFO,hdr_id);
+            }
+            else
+                rename_annots_push(args,ptr+2,str.s);
             *ptr = ':';
         }
         else if ( !strcasecmp("FILTER",str.s) )
@@ -2487,6 +2535,13 @@ static void init_columns(args_t *args)
                           "       (the annotation type is modified to \"Number=.\" and allele ordering is disregarded)\n");
                 fprintf(stderr,"Warning: the =INFO/TAG feature modifies the annotation to \"Number=.\" and disregards allele ordering\n");
             }
+
+            args->ncols++; args->cols = (annot_col_t*) realloc(args->cols,sizeof(annot_col_t)*args->ncols);
+            annot_col_t *col = &args->cols[args->ncols-1];
+            memset(col,0,sizeof(*col));
+            col->icol = icol;
+            col->replace = replace;
+
             int explicit_src_info = 0;
             int explicit_dst_info = 0;
             char *key_dst;
@@ -2517,15 +2572,14 @@ static void init_columns(args_t *args)
                     key_src[-2] = ':';
                     error("Did you mean \"FMT/%s\" rather than \"%s\"?\n",str.s,str.s);
                 }
+                else if ( !strcasecmp("FILTER",key_src) && args->tgts_is_vcf )
+                {
+                    col->getter = vcf_getter_filter2str;
+                }
             }
             else
                 key_src = key_dst;
 
-            args->ncols++; args->cols = (annot_col_t*) realloc(args->cols,sizeof(annot_col_t)*args->ncols);
-            annot_col_t *col = &args->cols[args->ncols-1];
-            memset(col,0,sizeof(*col));
-            col->icol = icol;
-            col->replace = replace;
             col->hdr_key_src = strdup(key_src);
             col->hdr_key_dst = strdup(key_dst);
 
@@ -2723,6 +2777,7 @@ static void rename_chrs(args_t *args, char *fname)
 // Dirty: this relies on bcf_hdr_sync NOT being called
 static int rename_annots_core(args_t *args, char *ori_tag, char *new_tag)
 {
+fprintf(stderr,"rename_annots_core ori=[%s] new=[%s]\n",ori_tag,new_tag);
     int type;
     if ( !strncasecmp("info/",ori_tag,5) ) type = BCF_HL_INFO, ori_tag += 5;
     else if ( !strncasecmp("format/",ori_tag,7) ) type = BCF_HL_FMT, ori_tag += 7;
@@ -2782,7 +2837,7 @@ static void rename_annots(args_t *args)
         while ( *ptr && isspace(*ptr) ) ptr++;
         if ( !*ptr ) { *rmme = ' '; error("Could not parse: %s\n", args->rename_annots_map[i]); }
         if ( rename_annots_core(args, args->rename_annots_map[i], ptr) < 0 )
-            error("Could not parse \"%s %s\", expected INFO, FORMAT, or FILTER prefix\n",args->rename_annots_map[i],ptr);
+            error("Cannot rename \"%s\" to \"%s\"\n",args->rename_annots_map[i],ptr);
     }
 }
 static void rename_annots_push(args_t *args, char *src, char *dst)
@@ -3084,6 +3139,8 @@ static int strstr_match(char *a, char *b)
 }
 static void annotate(args_t *args, bcf1_t *line)
 {
+    args->current_rec = line;
+
     int i, j;
     for (i=0; i<args->nrm; i++)
         args->rm[i].handler(args, line, &args->rm[i]);
@@ -3286,6 +3343,15 @@ static void annotate(args_t *args, bcf1_t *line)
             has_overlap = 1;
         }
     }
+    else if ( args->ncols )
+    {
+        for (j=0; j<args->ncols; j++)
+        {
+            if ( !args->cols[j].setter ) continue;
+            if ( args->cols[j].setter(args,line,&args->cols[j],NULL) )
+                error("fixme: Could not set %s at %s:%"PRId64"\n", args->cols[j].hdr_key_src,bcf_seqname(args->hdr,line),(int64_t) line->pos+1);
+        }
+    }
     if ( args->set_ids )
     {
         args->tmpks.l = 0;
