Extend per-population filtering

- Support per-population annotation for embedded functions, such as in
  -t 'MAC:1=int(MAF*AN)' -S groups.txt

- Fix custom functions to use rounding instead of truncation when converting to
  integer, so that -t 'MAC:1=float(MAF*AN)' and -t 'MAC:1=int(MAF*AN)' return
  consistent values.

Resolves #2514

Author: pd3

NEWS

@@ -23,6 +23,15 @@ Changes affecting specific commands:
     - Make --gvcf2vcf work around malformed gVCF records where INFO/END overlaps
       the next record
 
+* bcftools +fill-tags
+
+    - Fix custom functions to use rounding instead of truncation when converting to
+      integer, so that -t 'MAC:1=float(MAF*AN)' and -t 'MAC:1=int(MAF*AN)' return
+      consistent values.
+
+    - Support per-population annotation for embedded functions, such as in
+      -t 'MAC:1=int(MAF*AN)' -S groups.txt (#2514)
+
 * bcftools merge
 
     - Fixes in memory access with localized alleles, it could result in printing

filter.c

@@ -1,6 +1,6 @@
 /*  filter.c -- filter expressions.
 
-    Copyright (C) 2013-2025 Genome Research Ltd.
+    Copyright (C) 2013-2026 Genome Research Ltd.
 
     Author: Petr Danecek <pd3@sanger.ac.uk>
 
@@ -109,11 +109,22 @@ struct _filter_t
     float   *tmpf;
     kstring_t tmps;
     int max_unpack, mtmpi, mtmpf, nsamples;
+    uint8_t *sample_mask;   // bitmask restricting samples for FMT->INFO calculations
+    int sample_mask_set;    // number of samples after masking, 0: sample mask not used
     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;
+    struct {
+        bcf1_t *line;
+        int n_missing;
+        int n_samples;  // to be able to calculate F_MISSING
+        int AN;
+        int m_als;      // the size of allocated AC and AF
+        int32_t *AC;    // AC,AF are of size n_allele+1, idx=0 is for the REF allele, i.e. Number=R, not VCF's Number=A
+        double *AF;
+    } cached_site;
 #if ENABLE_PERL_FILTERS
     PerlInterpreter *perl;
 #endif
@@ -1443,26 +1454,81 @@ static void filters_set_nalt(filter_t *flt, bcf1_t *line, token_t *tok)
     tok->nvalues = 1;
     tok->values[0] = line->n_allele - 1;
 }
-static void filters_set_ac(filter_t *flt, bcf1_t *line, token_t *tok)
+#define CACHED flt->cached_site
+static int filters_cache_AC_stats(filter_t *flt, bcf1_t *line)
 {
-    hts_expand(int32_t, line->n_allele, flt->mtmpi, flt->tmpi);
-    if ( !bcf_calc_ac(flt->hdr, line, flt->tmpi, BCF_UN_INFO|BCF_UN_FMT) )
+    if ( CACHED.line==line ) return CACHED.n_samples > 0 ? 0 : -1;
+    CACHED.line = line;
+
+    if ( CACHED.m_als < line->n_allele )
     {
-        tok->nvalues = 0;
-        return;
+        CACHED.m_als = line->n_allele;
+        CACHED.AC = realloc(CACHED.AC,sizeof(*CACHED.AC)*CACHED.m_als);
+        CACHED.AF = realloc(CACHED.AF,sizeof(*CACHED.AF)*CACHED.m_als);
+    }
+
+    int i,j;
+    if ( !flt->sample_mask_set && bcf_calc_ac(flt->hdr, line, CACHED.AC, BCF_UN_INFO)>0 )
+    {
+        // all samples, can reuse INFO/AC,AN etc, when available
+        CACHED.n_samples = line->n_sample;
+        CACHED.n_missing = -1;  // cannot determine n_missing from INFO
+        CACHED.AN = 0;
+        for (i=0; i<line->n_allele; i++) CACHED.AN += CACHED.AC[i];
+        for (i=0; i<line->n_allele; i++) CACHED.AF[i] = CACHED.AN ? (double)CACHED.AC[i]/CACHED.AN : 0;
+        return 0;
+    }
+
+    CACHED.n_samples = flt->sample_mask_set ? flt->sample_mask_set : line->n_sample;
+    CACHED.n_missing = 0;
+    CACHED.AN = 0;
+    for (i=0; i<line->n_allele; i++) CACHED.AC[i] = 0;
+
+    int ngt = bcf_get_genotypes(flt->hdr, line, &flt->tmpi, &flt->mtmpi);
+    int ngt1 = ngt/line->n_sample;
+    for (i=0; i<line->n_sample; i++)
+    {
+        if ( flt->sample_mask_set && !flt->sample_mask[i] ) continue;
+        int32_t *ptr = flt->tmpi + ngt1*i;
+        int is_missing = 1;
+        for (j=0; j<ngt1; j++)
+        {
+            if ( bcf_gt_is_missing(ptr[j]) ) continue;
+            if ( ptr[j]==bcf_int32_vector_end ) break;
+            int ial = bcf_gt_allele(ptr[j]);
+            if ( ial >= line->n_allele )
+            {
+                static int warned = 0;
+                if ( !warned )
+                {
+                    fprintf(stderr,"The allele index too large (%d), skipping GT parsing at this site %s:%"PRId64". "
+                                   "(This warning is printed only once.)\n", ial,bcf_seqname(flt->hdr,line),line->pos+1);
+                    warned = 1;
+                }
+                CACHED.n_samples = 0;
+                return -1;
+            }
+            CACHED.AC[ial]++;
+            CACHED.AN++;
+            is_missing = 0;
+        }
+        if ( is_missing ) CACHED.n_missing++;
     }
-    int i, an = flt->tmpi[0];
-    for (i=1; i<line->n_allele; i++) an += flt->tmpi[i];
-    if ( !an )
+    for (i=0; i<line->n_allele; i++) CACHED.AF[i] = CACHED.AN ? (double)CACHED.AC[i]/CACHED.AN : 0;
+    return 0;
+}
+static void filters_set_ac(filter_t *flt, bcf1_t *line, token_t *tok)
+{
+    if ( filters_cache_AC_stats(flt,line) || !CACHED.AN )
     {
         tok->nvalues = 0;
         return;
     }
-    flt->tmpi[0] = an;  // for filters_set_[mac|af|maf]
+    // Here tok->idx is relative to ALT, i.e. interpreted as VCF's Number=A
     if ( tok->idx>=0 )
     {
         tok->nvalues = 1;
-        tok->values[0] = tok->idx+1<line->n_allele ? flt->tmpi[tok->idx+1] : 0;
+        tok->values[0] = tok->idx+1<line->n_allele ? CACHED.AC[tok->idx+1] : 0;
     }
     else if ( line->n_allele==1 )   // no ALT
     {
@@ -1471,45 +1537,68 @@ static void filters_set_ac(filter_t *flt, bcf1_t *line, token_t *tok)
     }
     else
     {
+        int i;
         hts_expand(double,line->n_allele,tok->mvalues,tok->values);
         for (i=1; i<line->n_allele; i++)
-            tok->values[i-1] = flt->tmpi[i];
+            tok->values[i-1] = CACHED.AC[i];
         tok->nvalues = line->n_allele - 1;
     }
 }
 static void filters_set_an(filter_t *flt, bcf1_t *line, token_t *tok)
 {
-    filters_set_ac(flt,line,tok);
-    tok->values[0] = tok->nvalues ? flt->tmpi[0] : 0;
     tok->nvalues = 1;
+    if ( filters_cache_AC_stats(flt,line)==0 )
+        tok->values[0] = CACHED.AN;
+    else
+        tok->values[0] = 0;
 }
 static void filters_set_mac(filter_t *flt, bcf1_t *line, token_t *tok)
 {
+    // Note: this one does not make much sense if nvalues!=1
     filters_set_ac(flt,line,tok);
     if ( !tok->nvalues ) return;
-    int i, an = flt->tmpi[0];
+    int i;
     for (i=0; i<tok->nvalues; i++)
-        if ( tok->values[i] > an*0.5 ) tok->values[i] = an - tok->values[i];
+        if ( tok->values[i] > 0.5*CACHED.AN ) tok->values[i] = CACHED.AN - tok->values[i];
 }
 static void filters_set_af(filter_t *flt, bcf1_t *line, token_t *tok)
 {
-    filters_set_ac(flt,line,tok);
-    if ( !tok->nvalues ) return;
-    int i, an = flt->tmpi[0];
-    for (i=0; i<tok->nvalues; i++)
-        tok->values[i] /= (double)an;
+    if ( filters_cache_AC_stats(flt,line) || !CACHED.AN )
+    {
+        tok->nvalues = 0;
+        return;
+    }
+    // Here tok->idx is relative to ALT, i.e. interpreted as VCF's Number=A
+    if ( tok->idx>=0 )
+    {
+        tok->nvalues = 1;
+        tok->values[0] = tok->idx+1<line->n_allele ? CACHED.AF[tok->idx+1] : 0;
+    }
+    else if ( line->n_allele==1 )   // no ALT
+    {
+        tok->nvalues = 1;
+        tok->values[0] = 0;
+    }
+    else
+    {
+        int i;
+        hts_expand(double,line->n_allele,tok->mvalues,tok->values);
+        for (i=1; i<line->n_allele; i++)
+            tok->values[i-1] = CACHED.AF[i];
+        tok->nvalues = line->n_allele - 1;
+    }
 }
 static void filters_set_maf(filter_t *flt, bcf1_t *line, token_t *tok)
 {
-    filters_set_ac(flt,line,tok);
+    filters_set_af(flt,line,tok);
     if ( !tok->nvalues ) return;
-    int i, an = flt->tmpi[0];
+    int i;
     for (i=0; i<tok->nvalues; i++)
     {
-        tok->values[i] /= (double)an;
         if ( tok->values[i] > 0.5 ) tok->values[i] = 1 - tok->values[i];
     }
 }
+#undef CACHED
 
 static int func_max(filter_t *flt, bcf1_t *line, token_t *rtok, token_t **stack, int nstack)
 {
@@ -4117,6 +4206,11 @@ static filter_t *filter_init_(bcf_hdr_t *hdr, const char *str, int exit_on_error
         }
     }
     filter->nsamples = filter->max_unpack&BCF_UN_FMT ? bcf_hdr_nsamples(filter->hdr) : 0;
+    if ( filter->nsamples )
+    {
+        filter->sample_mask = (uint8_t*) malloc(filter->nsamples);
+        if ( !filter->sample_mask ) error("Could not allocate %d bytes\n",filter->nsamples);
+    }
     for (i=0; i<nout; i++)
     {
         if ( out[i].tok_type==TOK_MAX )      { out[i].func = func_max; out[i].tok_type = TOK_FUNC; }
@@ -4193,12 +4287,15 @@ void filter_destroy(filter_t *filter)
     free(filter->used_tag);
     free(filter->cached_GT.buf);
     free(filter->cached_GT.mask);
+    free(filter->cached_site.AC);
+    free(filter->cached_site.AF);
     free(filter->filters);
     free(filter->flt_stack);
     free(filter->str);
     free(filter->tmpi);
     free(filter->tmpf);
     free(filter->tmps.s);
+    free(filter->sample_mask);
     free(filter);
 }
 
@@ -4389,12 +4486,20 @@ const double *filter_get_doubles(filter_t *filter, int *nval, int *nval1)
 
 void filter_set_samples(filter_t *filter, const uint8_t *samples)
 {
+    filter->sample_mask_set = 0;
+    if ( !filter->nsamples ) return;
+
     int i,j;
     for (i=0; i<filter->nfilters; i++)
     {
         if ( !filter->filters[i].nsamples ) continue;
         for (j=0; j<filter->filters[i].nsamples; j++) filter->filters[i].usmpl[j] = samples[j];
     }
+    for (i=0; i<filter->nsamples; i++)
+    {
+        filter->sample_mask[i] = samples[i];
+        if ( samples[i] ) filter->sample_mask_set++;
+    }
 }
 
 int filter_status(filter_t *filter)

plugins/fill-tags.c

@@ -1,6 +1,6 @@
 /* The MIT License
 
-   Copyright (c) 2015-2025 Genome Research Ltd.
+   Copyright (c) 2015-2026 Genome Research Ltd.
 
    Author: Petr Danecek <pd3@sanger.ac.uk>
 
@@ -283,7 +283,7 @@ void ftf_destroy(pop_t *pop)
 static inline int32_t int32_from_double(double src)
 {
     if (bcf_double_is_missing_or_vector_end(src)) return bcf_int32_missing;
-    else return (int)src;
+    else return (int)round(src);
 }
 
 int ftf_filter_expr(args_t *args, bcf1_t *rec, pop_t *pop, ftf_t *ftf)

test/fill-tags.MAC-MAF.1.out

@@ -0,0 +1,15 @@
+11	2343543	.	A	.	.	.	MAC_AB=0;MAC_BC=0;MAC=0	GT	0/0	0/0	0/0
+11	5464562	.	C	T	.	.	MAC_AB=.;MAC_BC=.;MAC=.	GT	./.	./.	./.
+20	76962	.	T	C	.	.	MAC_AB=1;MAC_BC=0;MAC=1	GT	0/1	1/1	1/1
+20	126310	.	ACC	A	.	.	MAC_AB=2;MAC_BC=1;MAC=2	GT	0/1	0/1	1/1
+20	138125	.	G	T	.	.	MAC_AB=2;MAC_BC=1;MAC=2	GT	0/1	0/1	1/1
+20	138148	.	C	T	.	.	MAC_AB=2;MAC_BC=1;MAC=2	GT	0/1	0/1	1/1
+20	271225	.	T	TTTA,TA	.	.	MAC_AB=1;MAC_BC=2;MAC=2	GT	0/2	0/1	1/2
+20	304568	.	C	T	.	.	MAC_AB=2;MAC_BC=1;MAC=2	GT	0|1	0|1	1|1
+20	326891	.	A	AC	.	.	MAC_AB=2;MAC_BC=1;MAC=2	GT	0|1	0|1	./.
+X	2928329	.	C	T	.	.	MAC_AB=0;MAC_BC=1;MAC=1	GT	0	0	0/1
+X	2933066	.	G	C	.	.	MAC_AB=0;MAC_BC=1;MAC=1	GT	0	0	0/1
+X	2942109	.	T	C	.	.	MAC_AB=1;MAC_BC=0;MAC=1	GT	0	1	1/1
+X	3048719	.	T	C	.	.	MAC_AB=1;MAC_BC=1;MAC=2	GT	0	1	0|1
+Y	8657215	.	C	A	.	.	MAC_AB=1;MAC_BC=0;MAC=1	GT	0	1	.
+Y	10011673	.	G	A	.	.	MAC_AB=0;MAC_BC=0;MAC=0	GT	1	1	.

test/fill-tags.MAC-MAF.vcf

@@ -0,0 +1,23 @@
+##fileformat=VCFv4.2
+##reference=file:///ref.fa
+##contig=<ID=11,length=135006516>
+##contig=<ID=20,length=63025520>
+##contig=<ID=X,length=155270560>
+##contig=<ID=Y,length=59373566>
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	NA00001	NA00002	NA00003
+11	2343543	.	A	.	.	.	.	GT	0/0	0/0	0/0
+11	5464562	.	C	T	.	.	.	GT	./.	./.	./.
+20	76962	.	T	C	.	.	.	GT	0/1	1/1	1/1
+20	126310	.	ACC	A	.	.	.	GT	0/1	0/1	1/1
+20	138125	.	G	T	.	.	.	GT	0/1	0/1	1/1
+20	138148	.	C	T	.	.	.	GT	0/1	0/1	1/1
+20	271225	.	T	TTTA,TA	.	.	.	GT	0/2	0/1	1/2
+20	304568	.	C	T	.	.	.	GT	0|1	0|1	1|1
+20	326891	.	A	AC	.	.	.	GT	0|1	0|1	./.
+X	2928329	.	C	T	.	.	.	GT	0	0	0/1
+X	2933066	.	G	C	.	.	.	GT	0	0	0/1
+X	2942109	.	T	C	.	.	.	GT	0	1	1/1
+X	3048719	.	T	C	.	.	.	GT	0	1	0|1
+Y	8657215	.	C	A	.	.	.	GT	0	1	.
+Y	10011673	.	G	A	.	.	.	GT	1	1	.

test/test.pl

@@ -694,6 +694,8 @@
 run_test(\&test_vcf_plugin,$opts,in=>'query.variantkey',out=>'query.add-variantkey.vcf',cmd=>'+add-variantkey',args=>'');
 run_test(\&test_vcf_plugin,$opts,in=>'query.variantkey',out=>'variantkey-hex.out',cmd=>'+variantkey-hex',args=>'test/');
 run_test(\&test_vcf_plugin,$opts,in=>'query.nucleotide',out=>'query.allele-length.tsv',cmd=>'+allele-length',args=>'');
+run_test(\&test_vcf_plugin,$opts,in=>'fill-tags.MAC-MAF',out=>'fill-tags.MAC-MAF.1.out',cmd=>'+fill-tags --no-version',args=>q[-- -S {PATH}/fill-tags.3.smpl -t 'MAC:1=MAC' | grep -v ^#]);
+run_test(\&test_vcf_plugin,$opts,in=>'fill-tags.MAC-MAF',out=>'fill-tags.MAC-MAF.1.out',cmd=>'+fill-tags --no-version',args=>q[-- -S {PATH}/fill-tags.3.smpl -t 'MAC:1=int(MAF*AN)' | grep -v ^#]);
 run_test(\&test_vcf_plugin,$opts,in=>'query.func.1',out=>'fill-tags.func.1.out',cmd=>'+fill-tags --no-version',args=>q[-- -t 'FMT/AB:1=float(FMT/AD[*:0] / ssum(FMT/AD[*]))' | grep -v ^#]);
 run_test(\&test_vcf_plugin,$opts,in=>'fmissing',out=>'fmissing.1.out',cmd=>'+fill-tags --no-version',args=>q[-- -S {PATH}/fmissing.txt -t 'F_MISSING' | grep -v ^#]);
 run_test(\&test_vcf_plugin,$opts,in=>'fmissing',out=>'fmissing.1.out',cmd=>'+fill-tags --no-version',args=>q[-- -S {PATH}/fmissing.txt -t 'F_MISSING:1=F_PASS(GT="mis")' | grep -v ^#]);