fix: multi-allelic statistics and capture filtering in variant-info

docs/advanced/coverage-evidence.md

@@ -17,12 +17,16 @@ into `N_HOM_REF`. The flags below decide *which* samples land there.
 ## What `N_NO_COVERAGE` represents
 
 `N_NO_COVERAGE` counts eligible samples whose hom-ref status is not trusted
-under the active criteria. The genotype invariant becomes:
+under the active criteria. At a biallelic site the genotype invariant becomes:
 
 ```
 N_HET + N_HOM_ALT + N_HOM_REF + N_FAIL + N_NO_COVERAGE = n_eligible
 ```
 
+At a multi-allelic site, eligible samples that carry only another ALT allele at
+the position fall in none of these categories for this allele, so the sum is
+lower than `n_eligible` by exactly that number of samples.
+
 Samples in `N_NO_COVERAGE` remain in `eligible` and contribute to `AN` (just
 like `N_FAIL`), so AC/AN/AF stay conservative — the field never inflates
 allele frequencies. Two rules always hold:
@@ -45,7 +49,7 @@ position. They run at query time, so no database rebuild is needed.
 
 | Flag | Effect |
 |------|--------|
-| `--min-pass K`     | A partially-covered tech must have ≥K PASS carriers (`het ∪ hom`) at the position. If it falls short, all of its non-carrier samples move from `N_HOM_REF` to `N_NO_COVERAGE`. |
+| `--min-pass K`     | A partially-covered tech must have ≥K PASS carriers (`het ∪ hom`) at the position. Carriers of any ALT allele there count, since each shows the position was sequenced. If it falls short, all of its non-carrier samples move from `N_HOM_REF` to `N_NO_COVERAGE`. |
 | `--min-observed K` | Same shape, but counts every recorded carrier (`het ∪ hom ∪ fail`). Useful when a non-PASS call still proves the position was sequenced. |
 
 When both flags are >0, both must hold (AND). The default `0` disables the
@@ -98,7 +102,7 @@ coverage decision is baked in.
 | `--min-dp D`     | Minimum `FORMAT/DP` per carrier. |
 | `--min-gq G`     | Minimum `FORMAT/GQ` per carrier. |
 | `--min-qual Q`   | Minimum VCF `QUAL` per carrier. |
-| `--min-covered K`| Per partially-covered tech, the position is "trusted" only if at least K of its carriers pass the quality thresholds. Non-carriers of failing positions are recorded as `N_NO_COVERAGE`. |
+| `--min-covered K`| Per partially-covered tech, the position is "trusted" only if at least K of its carriers pass the quality thresholds. Non-carriers of failing positions are recorded as `N_NO_COVERAGE`. The gate is evaluated at build time for each ALT allele separately; samples carrying another allele at the position are never counted as `N_NO_COVERAGE`. |
 
 A carrier counts as quality-passing only if **all** active thresholds hold
 (unset thresholds are simply ignored). At least one of these flags must be
@@ -133,7 +137,7 @@ afquery query --db ./db/ --locus chr1:925952 --min-quality-evidence 5
 ```
 
 `--min-quality-evidence K` requires each partially-covered tech to have ≥K
-quality-passing carriers at the position. Non-carriers of failing techs
+quality-passing carriers of any ALT allele at the position. Non-carriers of failing techs
 (other than those already filtered at build time) move to `N_NO_COVERAGE`.
 
 Running the flag against a database that was not built with quality data

docs/advanced/debugging-results.md

@@ -16,7 +16,7 @@ AN=0 means no eligible samples at the queried position. Work through these check
 | Position exists in database | `afquery query --db ./db/ --locus chr1:12345678` | If no result at all, the variant was not observed in any sample during ingestion. |
 | BED coverage (WES) | `afquery info --db ./db/` | If all eligible samples are WES and the position is outside capture regions, AN=0 is correct. |
 | Sample filter too restrictive | Remove `--phenotype` and `--sex` filters | Query with no filters first. If AN>0 without filters, the filter is excluding all samples. |
-| WES samples missing from AN | Compare `variant-info` carriers against `query` counts | If `variant-info` shows many WES carriers but AN reflects only the WGS samples, that technology's capture index is not matching the position. `afquery` warns on open when a capture BED is empty or matches no known chromosome, naming the offending contigs. |
+| WES samples missing from AN | `afquery query --db ./db/ --locus chr1:12345678 --tech WES_kit_A` | If this prints `No variants found for the given filters.` for a technology whose kit should cover the position, while the unfiltered query returns the variant, that technology's capture index is not matching the position. `afquery` warns on open when a capture BED is empty or matches no known chromosome, naming the offending contigs. |
 | AN=0 only on chrM | `afquery query --db ./db/ --locus chrM:3243` | Mitochondrial coverage is easy to lose on its own: a BED naming it `chrMT` still matches on the autosomes, so nothing else looks wrong. Names are normalized (`MT`, `M`, `chrMT` → `chrM`) since v0.4.0 — databases built before that may need the capture BED re-checked. |
 | Technology filter | Remove `--tech` filter | Check if any samples match the requested technology. |
 

docs/advanced/ploidy-and-sex-chroms.md

@@ -103,11 +103,11 @@ afquery query --db ./db/ --locus chrM:3243
 
 ### Counting Identity
 
-For every query result, the following identity holds:
+For every query result at a biallelic site, the following identity holds:
 
 **N_HET + N_HOM_ALT + N_HOM_REF + N_FAIL + N_NO_COVERAGE = n_eligible**
 
-This can be used to validate results. N_HOM_REF is the number of eligible samples that are homozygous reference (i.e., do not carry the alt allele and passed quality filters). N_NO_COVERAGE is 0 unless a coverage-evidence filter is active — see [Coverage Evidence](coverage-evidence.md).
+This can be used to validate results. At a multi-allelic site, eligible samples that carry only another ALT allele at the position fall in none of these categories, so the sum is lower by exactly that number of samples (see [Multi-allelic sites](../getting-started/understanding-output.md#multi-allelic-sites)). N_HOM_REF is the number of eligible samples that are homozygous reference (i.e., carry no ALT allele at the position and passed quality filters). N_NO_COVERAGE is 0 unless a coverage-evidence filter is active — see [Coverage Evidence](coverage-evidence.md).
 
 !!! note "Mutual exclusivity"
     N_HET, N_HOM_ALT, N_HOM_REF, N_FAIL, and N_NO_COVERAGE are mutually exclusive. A sample with a non-ref allele but FILTER≠PASS is counted in N_FAIL only — it does not appear in N_HET or N_HOM_ALT. Likewise, N_HOM_REF counts only PASS-filtered samples.

docs/getting-started/preprocessing.md

@@ -2,7 +2,15 @@
 
 AFQuery ingests single-sample, normalized VCF files. While AFQuery itself does not perform VCF normalization, the accuracy of your allele frequency estimates depends on the quality and consistency of your input VCFs. This page explains why normalization matters and provides a reference pipeline.
 
-A ready-to-use normalization script is provided at [`resources/normalize_vcf.sh`](https://github.com/babelomics/afquery/blob/master/resources/normalize_vcf.sh).
+A ready-to-use normalization script is provided at [`resources/normalize_vcf.sh`](https://github.com/babelomics/afquery/blob/master/resources/normalize_vcf.sh). It requires **bcftools 1.20 or later**: older releases refuse to split multi-allelic records and remove duplicates in the same `bcftools norm` call, and the script exits with an error if it finds one.
+
+!!! warning "VCFs normalized with the script from afquery 0.4.2 or earlier"
+    Those versions of the script did not split multi-allelic records, and removed
+    one of two different indels at the same position as a duplicate. Indels from
+    multi-allelic sites could therefore be stored untrimmed (e.g. `CAAAAAAA>CAAAAAAAAAAAA`
+    instead of `C>CAAAAA`), so the same variant appeared under two representations, or
+    were dropped. Re-normalize those VCFs with the current script and rebuild the
+    database to correct them.
 
 
 ---
@@ -17,7 +25,7 @@ Multi-allelic variants and complex indels can be represented in multiple equival
 - Multi-allelic sites may not be decomposed into biallelic records
 - Duplicate records can inflate AC
 
-`bcftools norm` left-aligns indels against the reference genome and decomposes multi-allelic sites, ensuring consistent representation across samples.
+`bcftools norm -m -both` decomposes multi-allelic sites into biallelic records, and left-aligns and trims each allele against the reference genome, ensuring consistent representation across samples. `-d exact` then drops only records whose alleles are identical, so two different alleles at the same position are both kept.
 
 ### Ploidy correction for sex chromosomes
 

docs/getting-started/understanding-output.md

@@ -13,7 +13,7 @@ This page explains what each field in AFQuery output means and how to interpret
 | **AF** | float | Allele frequency — `AC / AN`. `None` when AN=0 |
 | **N_HET** | int | Number of eligible samples heterozygous for the alt allele (GT=0/1) |
 | **N_HOM_ALT** | int | Number of eligible samples homozygous for the alt allele (GT=1/1 or GT=1). Includes haploid carriers on sex chromosomes and chrM. See [Ploidy](../advanced/ploidy-and-sex-chroms.md#genotype-counting). |
-| **N_HOM_REF** | int | Number of eligible samples homozygous reference (GT=0/0 or GT=0) |
+| **N_HOM_REF** | int | Number of eligible samples homozygous reference (GT=0/0 or GT=0). Samples carrying another ALT allele at the same position are not counted (see [Multi-allelic sites](#multi-allelic-sites)). |
 | **n_eligible** | int | Number of eligible samples — those passing the sex/phenotype/tech filters *and* covered at this position |
 | **N_FAIL** | int | Number of eligible samples whose call at this position had FILTER≠PASS. These samples are counted *only* in N_FAIL — not in N_HET, N_HOM_ALT, or N_HOM_REF — but they stay eligible and still count toward AN. |
 | **N_NO_COVERAGE** | int | Number of eligible samples whose tech lacks coverage evidence at this position. Excluded from `N_HOM_REF` to keep AC/AN conservative. Always `0` unless a coverage-evidence filter is active. See [Coverage Evidence](../advanced/coverage-evidence.md). |
@@ -84,6 +84,19 @@ AN=0 means no eligible samples have coverage at this position. This happens when
 !!! warning "AN=0 does not mean the variant is absent"
     AN=0 means AFQuery has no data to compute frequency. It is not evidence of rarity.
 
+### Multi-allelic sites
+
+Each ALT allele at a position is reported on its own line, and its counts refer to
+that allele only. A sample carrying a different ALT allele at the same position
+(for example `G/T` on the `G>A` line) is neither a carrier nor homozygous
+reference for this allele, so it appears in none of `N_HET`, `N_HOM_ALT`,
+`N_HOM_REF`, `N_FAIL` or `N_NO_COVERAGE`. The five counts add up to
+`n_eligible` at biallelic sites; at a multi-allelic site they fall short by the
+number of such samples. AC, AN and AF are not affected.
+
+A sample with two different ALT alleles (`GT=1/2`) is counted as heterozygous on
+both lines.
+
 ### Warnings
 
 afquery emits a `AfqueryWarning` to stderr when a query may silently return fewer or no results. Common causes:
@@ -126,7 +139,7 @@ When using `afquery annotate`, the following INFO fields are added to each varia
 | `AFQUERY_AF` | A (per ALT) | Allele frequency — one value per ALT allele |
 | `AFQUERY_N_HET` | A (per ALT) | Heterozygous sample count per ALT allele |
 | `AFQUERY_N_HOM_ALT` | A (per ALT) | Homozygous alt sample count per ALT allele |
-| `AFQUERY_N_HOM_REF` | A (per ALT) | Homozygous ref sample count per ALT allele |
+| `AFQUERY_N_HOM_REF` | A (per ALT) | Homozygous ref sample count per ALT allele, excluding carriers of other ALT alleles at the position |
 | `AFQUERY_N_FAIL` | 1 (per site) | Fail sample count — shared across all ALT alleles |
 | `AFQUERY_N_NO_COVERAGE` | A (per ALT) | Eligible samples whose tech lacks coverage evidence at this position. Always `0` unless a coverage-evidence filter is active. See [Coverage Evidence](../advanced/coverage-evidence.md). |
 

docs/guides/annotate-vcf.md

@@ -25,7 +25,7 @@ afquery annotate \
 | `AFQUERY_AF` | Float | A (per ALT) | Allele frequency (`AC / AN`) |
 | `AFQUERY_N_HET` | Integer | A (per ALT) | Heterozygous sample count |
 | `AFQUERY_N_HOM_ALT` | Integer | A (per ALT) | Homozygous alt sample count |
-| `AFQUERY_N_HOM_REF` | Integer | A (per ALT) | Homozygous ref sample count |
+| `AFQUERY_N_HOM_REF` | Integer | A (per ALT) | Homozygous ref sample count, excluding carriers of other ALT alleles at the position (including alleles stored in the database but absent from the input record) |
 | `AFQUERY_N_FAIL` | Integer | 1 (per site) | Eligible samples whose call had FILTER≠PASS. Excluded from AC, but still counted in AN. Mutually exclusive with N_HET/N_HOM_ALT/N_HOM_REF. |
 | `AFQUERY_N_NO_COVERAGE` | Integer | A (per ALT) | Eligible samples whose tech lacks coverage evidence at this position. Excluded from `N_HOM_REF` to keep AC/AN conservative. Always `0` unless a coverage-evidence filter is active. See [Coverage Evidence](../advanced/coverage-evidence.md). |
 

docs/guides/query.md

@@ -155,7 +155,7 @@ that fall below a threshold are reported in **N_NO_COVERAGE** instead of N_HOM_R
 
 | Flag | Meaning |
 |------|---------|
-| `--min-pass K` | A partially-covered tech is valid for hom-ref at a position only if it has ≥K PASS carriers (het\|hom). Otherwise its non-carrier samples move to `N_NO_COVERAGE`. |
+| `--min-pass K` | A partially-covered tech is valid for hom-ref at a position only if it has ≥K PASS carriers (het\|hom) of any ALT allele there. Otherwise its non-carrier samples move to `N_NO_COVERAGE`. |
 | `--min-observed K` | Same as `--min-pass`, but counts any VCF entry (`het\|hom\|fail`). Useful when you want to include calls that failed FILTER as evidence the position was sequenced. |
 | `--min-quality-evidence K` | Requires ≥K quality-passing carriers per partially-covered tech. Requires a database built with `--min-dp`, `--min-gq`, `--min-qual`, or `--min-covered`. |
 
@@ -167,8 +167,11 @@ afquery query --db ./db/ --locus chr1:925952 --min-pass 1
 afquery query --db ./db/ --region chr1:900000-1000000 --min-observed 2 --min-pass 1
 ```
 
-The genotype invariant becomes:
+At a biallelic site the genotype invariant becomes:
 `N_HET + N_HOM_ALT + N_HOM_REF + N_FAIL + N_NO_COVERAGE = n_eligible`.
+At a multi-allelic site, eligible samples that carry only another ALT allele at
+the position fall in none of these categories for this allele, so the sum is
+lower than `n_eligible` by exactly that number of samples.
 
 Fully-covered samples (those whose tech was registered without a BED) are
 never affected. Carrier samples (het/hom/fail) are never moved to

docs/guides/variant-info.md

@@ -15,17 +15,14 @@ afquery variant-info --db ./db/ --locus chr1:925952
 !!! tip
     `variant-info` is the natural next step after `query` — once you find a variant of interest, use it to see which specific samples carry it.
 
-!!! note "Carrier counts may exceed the counts reported by `query`"
-    `variant-info` lists every carrier present in the source VCFs. `query` additionally
-    restricts to samples whose capture regions cover the position, so a WES sample with a
-    call just outside its own BED appears here but is not counted in `AN`, `AC` or the
-    genotype tallies. The two commands answer different questions — "who carries it?"
-    versus "what is the frequency among samples that could have been called?" — so a
-    modest difference is expected.
-
-    A *large* gap is worth investigating: if `variant-info` shows many carriers from a
-    technology that contributes nothing to `AN`, check that technology's capture index
-    (see [Debugging Results](../advanced/debugging-results.md)).
+!!! note "Carriers match the counts reported by `query`"
+    `variant-info` lists only samples whose capture regions cover the position, the same
+    eligible set `query` uses. A WES sample with a call outside its own BED is not counted
+    in `AN`, `AC` or the genotype tallies, and is not listed here either. On diploid
+    regions the number of `het`, `hom` and `alt` carriers therefore equals `N_HET`,
+    `N_HOM_ALT` and `N_FAIL` for the same allele and filters. (In haploid regions a
+    single-allele call is listed as `het` but counted in `N_HOM_ALT`.) Earlier versions
+    also listed off-target carriers, so the two commands could disagree.
 
 By default all samples are queried and results are printed as an aligned text table:
 

docs/reference/python-api.md

@@ -457,8 +457,11 @@ class QueryResult:
     N_NO_COVERAGE: int       # Eligible samples whose tech lacks evidence (excluded from N_HOM_REF)
 ```
 
-The new genotype invariant is
+At a biallelic site the genotype categories partition the eligible samples:
 `N_HET + N_HOM_ALT + N_HOM_REF + N_FAIL + N_NO_COVERAGE == n_samples_eligible`.
+At a multi-allelic site, eligible samples that carry only another ALT allele at
+the position fall in none of these categories for this allele, so the sum is
+lower than `n_samples_eligible` by exactly that number of samples.
 See [Coverage Evidence](../advanced/coverage-evidence.md) for details on
 `N_NO_COVERAGE`.
 

resources/normalize_vcf.sh

@@ -25,6 +25,17 @@ REF="$3"
 GENDER="$4"
 THREADS="${5:-4}"
 
+# -----------------------------
+# bcftools version (< 1.20 refuses to combine -m and -d in bcftools norm)
+# -----------------------------
+MIN_BCFTOOLS="1.20"
+BCFTOOLS_VERSION=$(bcftools --version | head -n1 | awk '{print $2}')
+if [[ "$(printf '%s\n' "${MIN_BCFTOOLS}" "${BCFTOOLS_VERSION}" | sort -V | head -n1)" != "${MIN_BCFTOOLS}" ]]
+then
+    echo "ERROR: bcftools >= ${MIN_BCFTOOLS} is required (found ${BCFTOOLS_VERSION})"
+    exit 1
+fi
+
 OUT="${VCF_ID}_norm.vcf.gz"
 CHR_MAP="chr_list.txt"
 GENDER_FILE="${VCF_ID}_gender.txt"
@@ -88,8 +99,8 @@ bcftools annotate \
     --rename-chrs "${CHR_MAP}" \
     "${VCF}" | \
 bcftools norm \
+    -m -both \
     -d exact \
-    -d both \
     -f "${REF}" \
     --check-ref ws \
     --targets "${TARGETS}" | \