Arrays in Number=R tags can be now subscripted by alleles found in FORMAT/GT

For example

  FORMAT/AD[GT] > 10        .. require support of more than 10 reads for each allele
  FORMAT/AD[0:GT] > 10      .. same as above, but in the first sample
  sSUM(FORMAT/AD[GT]) > 20  .. require total sample depth bigger than 20

Resolves #1849

Author: pd3

NEWS

@@ -13,6 +13,12 @@ Changes affecting the whole of bcftools, or multiple commands:
         - evaluate to C_i = A_0 + B_i when length(A)=1 and set length(C)=length(B)
         - throw an error when length(A)!=length(B) AND length(A)!=1 AND length(B)!=1
 
+    - Arrays in Number=R tags can be now subscripted by alleles found in FORMAT/GT. For example,
+
+        FORMAT/AD[GT] > 10        .. require support of more than 10 reads for each allele
+        FORMAT/AD[0:GT] > 10      .. same as above, but in the first sample
+        sSUM(FORMAT/AD[GT]) > 20  .. require total sample depth bigger than 20
+
 
 Changes affecting specific commands:
 

doc/bcftools.txt

@@ -3505,7 +3505,7 @@ to require that all alleles are of the given type. Compare
 
 * array subscripts (0-based), "*" for any element, "-" to indicate a range. Note that
 for querying FORMAT vectors, the colon ":" can be used to select a sample and an
-element of the vector, as shown in the examples below
+element of the vector, as shown in the examples below.
 
         INFO/AF[0] > 0.3             .. first AF value bigger than 0.3
         FORMAT/AD[0:0] > 30          .. first AD value of the first sample bigger than 30
@@ -3522,6 +3522,14 @@ element of the vector, as shown in the examples below
         (DP4[0]+DP4[1])/(DP4[2]+DP4[3]) > 0.3
         CSQ[*] ~ "missense_variant.*deleterious"
 
+* in addition to array subscripts shown above, it is possible to subscript arrays
+of Number=R tags by alleles found in FORMAT/GT (starting with version 1.17).
+For example
+
+        FORMAT/AD[GT] > 10        .. require support of more than 10 reads for each allele
+        FORMAT/AD[0:GT] > 10      .. same as above, but in the first sample
+        sSUM(FORMAT/AD[GT]) > 20  .. require total sample depth bigger than 20
+
 * with many samples it can be more practical to provide a file with sample names,
 one sample name per line
 

filter.c

@@ -69,6 +69,7 @@ typedef struct _token_t
     int hdr_id, tag_type;   // BCF header lookup ID and one of BCF_HL_* types
     int idx;            // 0-based index to VCF vectors,
                         //  -2: list (e.g. [0,1,2] or [1..3] or [1..] or any field[*], which is equivalent to [0..])
+                        //  -3: select indices on the fly based on values in GT
     int *idxs;          // set indexes to 0 to exclude, to 1 to include, and last element negative if unlimited; used by VCF retrievers only
     int nidxs, nuidxs;  // size of idxs array and the number of elements set to 1
     uint8_t *usmpl;     // bitmask of used samples as set by idx, set for FORMAT fields, NULL otherwise
@@ -100,6 +101,11 @@ struct _filter_t
     float   *tmpf;
     kstring_t tmps;
     int max_unpack, mtmpi, mtmpf, nsamples;
+    struct {
+        bcf1_t *line;
+        int32_t *buf, nbuf, mbuf;   // GTs as obtained by bcf_get_genotypes()
+        uint64_t *mask;             // GTs as mask, e.g 0/0 is 1; 0/1 is 3, max 63 unique alleles
+    } cached_GT;
 #if ENABLE_PERL_FILTERS
     PerlInterpreter *perl;
 #endif
@@ -350,6 +356,44 @@ char *expand_path(char *path)
 
     return strdup(path);
 }
+static int filters_cache_genotypes(filter_t *flt, bcf1_t *line)
+{
+    if ( flt->cached_GT.line==line ) return flt->cached_GT.nbuf > 0 ? 0 : -1;
+    flt->cached_GT.line = line;
+    flt->cached_GT.nbuf = bcf_get_genotypes(flt->hdr, line, &flt->cached_GT.buf, &flt->cached_GT.mbuf);
+    if ( flt->cached_GT.nbuf<=0 ) return -1;
+    if ( !flt->cached_GT.mask )
+    {
+        flt->cached_GT.mask = (uint64_t*) malloc(sizeof(*flt->cached_GT.mask)*flt->nsamples);
+        if ( !flt->cached_GT.mask ) error("Could not alloc %zu bytes\n",sizeof(*flt->cached_GT.mask)*flt->nsamples);
+    }
+    int i,j, ngt1 = flt->cached_GT.nbuf / line->n_sample;
+    for (i=0; i<line->n_sample; i++)
+    {
+        int32_t *ptr = flt->cached_GT.buf + i*ngt1;
+        flt->cached_GT.mask[i] = 0;
+        for (j=0; j<ngt1; j++)
+        {
+            if ( bcf_gt_is_missing(ptr[j]) ) continue;
+            if ( ptr[j]==bcf_int32_vector_end ) break;
+            int allele = bcf_gt_allele(ptr[j]);
+            if ( allele > 63 )
+            {
+                static int warned = 0;
+                if ( !warned )
+                {
+                    fprintf(stderr,"Too many alleles, skipping GT filtering at this site %s:%"PRId64". "
+                                   "(This warning is printed only once.)\n", bcf_seqname(flt->hdr,line),line->pos+1);
+                    warned = 1;
+                }
+                flt->cached_GT.nbuf = 0;
+                return -1;
+            }
+            flt->cached_GT.mask[i] |= 1<<allele;
+        }
+    }
+    return 0;
+}
 
 static void filters_set_qual(filter_t *flt, bcf1_t *line, token_t *tok)
 {
@@ -762,6 +806,27 @@ static void filters_set_format_int(filter_t *flt, bcf1_t *line, token_t *tok)
                 tok->values[i] = ptr[tok->idx];
         }
     }
+    else if ( tok->idx==-3 )
+    {
+        if ( filters_cache_genotypes(flt,line)!=0 )
+        {
+            tok->nvalues = 0;
+            return;
+        }
+        for (i=0; i<tok->nsamples; i++)
+        {
+            if ( !tok->usmpl[i] ) continue;
+            int32_t *src = flt->tmpi + i*nsrc1;
+            double  *dst = tok->values + i*tok->nval1;
+            int k, j = 0;
+            for (k=0; k<nsrc1; k++) // source values are AD[0..nsrc1]
+            {
+                if ( !(flt->cached_GT.mask[i] & (1<<k)) ) continue;
+                dst[j++] = src[k];
+            }
+            for (; j<tok->nval1; j++) bcf_double_set_vector_end(dst[j]);
+        }
+    }
     else
     {
         int kend = tok->idxs[tok->nidxs-1] < 0 ? tok->nval1 : tok->nidxs;
@@ -825,6 +890,33 @@ static void filters_set_format_float(filter_t *flt, bcf1_t *line, token_t *tok)
                 tok->values[i] = ptr[tok->idx];
         }
     }
+    else if ( tok->idx==-3 )
+    {
+        if ( filters_cache_genotypes(flt,line)!=0 )
+        {
+            tok->nvalues = 0;
+            return;
+        }
+        for (i=0; i<tok->nsamples; i++)
+        {
+            if ( !tok->usmpl[i] ) continue;
+            float  *src = flt->tmpf + i*nsrc1;
+            double *dst = tok->values + i*tok->nval1;
+            int k, j = 0;
+            for (k=0; k<nsrc1; k++) // source values are AF[0..nsrc1]
+            {
+                if ( !(flt->cached_GT.mask[i] & (1<<k)) ) continue;
+                if ( bcf_float_is_missing(src[k]) )
+                    bcf_double_set_missing(dst[j]);
+                else if ( bcf_float_is_vector_end(src[k]) )
+                    bcf_double_set_vector_end(dst[j]);
+                else
+                    dst[j] = src[k];
+                j++;
+            }
+            for (; j<tok->nval1; j++) bcf_double_set_vector_end(dst[j]);
+        }
+    }
     else
     {
         int kend = tok->idxs[tok->nidxs-1] < 0 ? tok->nval1 : tok->nidxs;
@@ -989,9 +1081,9 @@ static void _filters_set_genotype(filter_t *flt, bcf1_t *line, token_t *tok, int
     tok->str_value.s[tok->str_value.l] = 0;
     tok->nval1 = nvals1;
 }
-static void filters_set_genotype2(filter_t *flt, bcf1_t *line, token_t *tok) { _filters_set_genotype(flt, line, tok, 2); }
-static void filters_set_genotype3(filter_t *flt, bcf1_t *line, token_t *tok) { _filters_set_genotype(flt, line, tok, 3); }
-static void filters_set_genotype4(filter_t *flt, bcf1_t *line, token_t *tok) { _filters_set_genotype(flt, line, tok, 4); }
+static void filters_set_genotype2(filter_t *flt, bcf1_t *line, token_t *tok) { _filters_set_genotype(flt, line, tok, 2); }  // rr, ra, aa, aA etc
+static void filters_set_genotype3(filter_t *flt, bcf1_t *line, token_t *tok) { _filters_set_genotype(flt, line, tok, 3); }  // hap, hom, het
+static void filters_set_genotype4(filter_t *flt, bcf1_t *line, token_t *tok) { _filters_set_genotype(flt, line, tok, 4); }  // mis, alt, ref
 
 static void filters_set_genotype_string(filter_t *flt, bcf1_t *line, token_t *tok)
 {
@@ -2480,6 +2572,14 @@ static int parse_idxs(char *tag_idx, int **idxs, int *nidxs, int *idx)
         *idx       = -2;
         return 0;
     }
+    if ( !strcmp("GT", tag_idx) )
+    {
+        *idxs = (int*) malloc(sizeof(int));
+        (*idxs)[0] = -1;
+        *nidxs     = 1;
+        *idx = -3;
+        return 0;
+    }
 
     // TAG[integer] .. one field; idx positive
     char *end, *beg = tag_idx;
@@ -2595,7 +2695,7 @@ static void parse_tag_idx(bcf_hdr_t *hdr, int is_fmt, char *tag, char *tag_idx,
                 tok->idxs = (int*) malloc(sizeof(int));
                 tok->idxs[0] = -1;
                 tok->nidxs   = 1;
-                tok->idx     = -2;
+                tok->idx     = idx1;
             }
             else if ( bcf_hdr_id2number(hdr,BCF_HL_FMT,tok->hdr_id)!=1 )
                 error("The FORMAT tag %s can have multiple subfields, run as %s[sample:subfield]\n", tag,tag);
@@ -2620,7 +2720,7 @@ static void parse_tag_idx(bcf_hdr_t *hdr, int is_fmt, char *tag, char *tag_idx,
             if ( idx1 >= bcf_hdr_nsamples(hdr) ) error("The sample index is too large: %s\n", ori);
             tok->usmpl[idx1] = 1;
         }
-        else if ( idx1==-2 )
+        else if ( idx1==-2 || idx1==-3 )
         {
             for (i=0; i<nidxs1; i++)
             {
@@ -2820,7 +2920,11 @@ static int filters_init1(filter_t *filter, char *str, int len, token_t *tok)
         if ( is_fmt==-1 ) is_fmt = 0;
     }
     if ( is_array )
+    {
         parse_tag_idx(filter->hdr, is_fmt, tmp.s, tmp.s+is_array, tok);
+        if ( tok->idx==-3 && bcf_hdr_id2length(filter->hdr,BCF_HL_FMT,tok->hdr_id)!=BCF_VL_R )
+            error("Error: GT subscripts can be used only with Number=R tags\n");
+    }
     else if ( is_fmt && !tok->nsamples )
     {
         int i;
@@ -3525,6 +3629,8 @@ void filter_destroy(filter_t *filter)
             free(filter->filters[i].regex);
         }
     }
+    free(filter->cached_GT.buf);
+    free(filter->cached_GT.mask);
     free(filter->filters);
     free(filter->flt_stack);
     free(filter->str);

test/query.92.out

@@ -0,0 +1,6 @@
+3000151	0/0	A	10,0,0
+3000151	0/0	B	10,0,0
+3000155	0/2	A	20,0,10
+3000155	0/2	B	200,0,10
+3000157	1/2	A	0,40,10
+3000157	1/2	B	0,400,10

test/query.93.out

@@ -0,0 +1,2 @@
+3000153	2/2	B	0,0,210
+3000155	0/2	B	200,0,10

test/query.94.out

@@ -0,0 +1,3 @@
+3000154	0/1	A	30,13,0
+3000156	1/0	A	30,14,0
+3000158	2/0	A	0,5,30

test/setGT.6.2.out

@@ -0,0 +1,7 @@
+##fileformat=VCFv4.2
+##FILTER=<ID=PASS,Description="All filters passed">
+##contig=<ID=1,assembly=b37,length=249250621>
+##reference=file:ref.fa
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	A	B	C	D
+1	1890	.	A	G	.	.	.	GT	./.	1|.	./.	1|1

test/test.pl

@@ -184,6 +184,9 @@
 run_test(\&test_vcf_query,$opts,in=>'query.filter.5',out=>'query.57.out',args=>q[-f'[%POS\\t%SAMPLE\\t%GT\\t%AD\\n]' -i'GT="het" & binom(FMT/AD)>0.01']);
 run_test(\&test_vcf_query,$opts,in=>'query.filter.5',out=>'query.58.out',args=>q[-f'[%POS\\t%SAMPLE\\t%GT\\t%AD\\n]' -i'GT="het" & binom(FMT/AD[:0],FMT/AD[:1])>0.01']);
 run_test(\&test_vcf_query,$opts,in=>'query.filter.5',out=>'query.59.out',args=>q[-f'%POS\\t%AD\\n' -i'binom(INFO/AD[0],INFO/AD[1])>0.01']);
+run_test(\&test_vcf_query,$opts,in=>'query.filter.5',out=>'query.92.out',args=>q[-f'[%POS\\t%GT\\t%SAMPLE\\t%AD\\n]' -i'FMT/AD[GT]==10']);
+run_test(\&test_vcf_query,$opts,in=>'query.filter.5',out=>'query.93.out',args=>q[-f'[%POS\\t%GT\\t%SAMPLE\\t%AD\\n]' -i'sSUM(FMT/AD[GT])==210']);
+run_test(\&test_vcf_query,$opts,in=>'query.filter.5',out=>'query.94.out',args=>q[-f'[%POS\\t%GT\\t%SAMPLE\\t%AD\\n]' -i'FMT/AD[0:GT]==30']);
 run_test(\&test_vcf_query,$opts,in=>'query',out=>'query.60.out',args=>q[-f'%CHROM %POS\\n' -i'CHROM="4"']);
 run_test(\&test_vcf_query,$opts,in=>'query.negative',out=>'query.61.out',args=>q[-f'%POS\\t%TAG1\\n' -i'(TAG1>=-129 && TAG1<=-120) || (TAG1>=-32769 && TAG1<=-32760)']);
 run_test(\&test_vcf_query,$opts,in=>'query.negative',out=>'query.61.out',args=>q[-f'%POS\\t%TAGV1\\n' -i'(TAGV1>=-129 && TAGV1<=-120) || (TAGV1>=-32769 && TAGV1<=-32760)']);