Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -158,7 +158,7 @@ Every query and annotated VCF reports:
| `N_HET` | Heterozygous PASS carriers |
| `N_HOM_ALT` | Homozygous-alt PASS carriers |
| `N_HOM_REF` | Samples trusted as homozygous reference |
| `N_FAIL` | Carriers with FILTER ≠ PASS (excluded from AC/AN) |
| `N_FAIL` | Samples with a failed call (FILTER ≠ PASS) at the position (excluded from AC/AN) |
| `N_NO_COVERAGE` | Non-carriers on a partially-covered tech without sufficient evidence to call hom-ref |

VCF INFO field names use the `AFQUERY_` prefix (e.g. `AFQUERY_AF`, `AFQUERY_N_NO_COVERAGE`).
Expand Down
12 changes: 8 additions & 4 deletions docs/advanced/coverage-evidence.md
Original file line number Diff line number Diff line change
Expand Up @@ -27,9 +27,8 @@ 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:
Samples in `N_NO_COVERAGE` remain in `eligible` and contribute to `AN`, so
AC/AN/AF stay conservative — the field never inflates allele frequencies. Two rules always hold:

- **Carriers are never reclassified.** A sample with a `het`, `hom`, or
`fail` call at the position stays in its category. `N_NO_COVERAGE` only
Expand All @@ -55,6 +54,11 @@ position. They run at query time, so no database rebuild is needed.
When both flags are >0, both must hold (AND). The default `0` disables the
gate.

`afquery annotate` applies the same gates to input alleles that are not in the
database. At a position where no sample has any call, no tech has evidence, so
under any active gate the non-carriers of partially-covered techs are reported
as `N_NO_COVERAGE`.

!!! tip
If your VCFs do not carry `FORMAT/DP` or `FORMAT/GQ`, these are the
flags you want. They are the cheapest option and apply to any database.
Expand Down Expand Up @@ -102,7 +106,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`. 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`. |
| `--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`. Carriers of any ALT allele at the position count, so every allele of a multi-allelic site gets the same verdict; samples carrying another allele there 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
Expand Down
4 changes: 2 additions & 2 deletions docs/advanced/debugging-results.md
Original file line number Diff line number Diff line change
Expand Up @@ -42,13 +42,13 @@ A variant you expect to find is not in the database:

### 4. Unexpected N_FAIL > 0

`N_FAIL > 0` means some eligible samples had a call with `FILTER≠PASS` at this position. These samples are excluded from AC, but they remain eligible and still count in AN. This is usually benign (1–2 samples), but a high N_FAIL warrants investigation:
`N_FAIL > 0` means some eligible samples had a call with `FILTER≠PASS` at this position. These samples are left out of AC and AN. This is usually benign (1–2 samples), but a high N_FAIL warrants investigation:

| N_FAIL relative to n_eligible | Likely cause |
|---|---|
| 1–2 samples | Isolated low-quality calls — not concerning |
| > 5% of n_eligible | Systematic sequencing artifact at this site |
| All eligible samples | Site-wide QC failure — AF=0 but variant is present in source VCFs |
| All eligible samples | Site-wide QC failure — AN=0 and AF is `NA`, although the variant is present in source VCFs |

To inspect a site with high N_FAIL, query with `--format json` to see all fields:

Expand Down
93 changes: 32 additions & 61 deletions docs/advanced/filter-pass-tracking.md
Original file line number Diff line number Diff line change
@@ -1,106 +1,77 @@
# FILTER=PASS Tracking

AFQuery tracks variants that are called but fail quality filters (`FILTER≠PASS`). This allows distinguishing between "variant not present" and "variant present but low quality" in your cohort.
AFQuery keeps the calls that fail quality filters (`FILTER≠PASS`) instead of discarding them. A failed call tells you the position was looked at but the sample's genotype there is unknown, which is different both from "variant not present" and from "variant present".

---

## Background: VCF FILTER Field
## What counts as a failed call

In VCF format, the FILTER column indicates whether a variant call passed quality filters:
A sample has a failed call at a position when its VCF holds, for any ALT allele there:

- `PASS` or `.` (missing) — the variant passed all filters
- Any other value (e.g., `LowQual`, `VQSRTrancheSNP99.90to100.00`) — the variant failed one or more filters
- a call carrying the ALT allele with `FILTER` other than `PASS` or `.`;
- a missing genotype (`./.`) at a site that failed a filter; or
- calls that carry more alleles than the sample has at the position, such as two ALT alleles on chrX in a male. These are detected when the VCF is ingested, which logs a warning (see [Genotypes inconsistent with ploidy](ploidy-and-sex-chroms.md#genotypes-inconsistent-with-ploidy)).

AFQuery's default behavior:

- Only `FILTER=PASS` calls contribute to **AC** (and therefore to **AF**). This is always enforced.
- **AN** is not affected by the filter — it counts every eligible sample at the position, failed calls included.
These samples are stored in the variant's `fail_bitmap` and reported as `N_FAIL`. A homozygous-reference call (`0/0`) is never stored, so at a failed site it would count as hom-ref; [`normalize_vcf.sh`](../getting-started/preprocessing.md) avoids this by masking every genotype at a failed site as missing.

---

## fail_bitmap

AFQuery stores a third bitmap per variant alongside `het_bitmap` and `hom_bitmap`:

- **`fail_bitmap`** — bit set for each sample whose call at this site has `FILTER≠PASS`. In practice this is two cases: a sample that carries the alt allele but whose call failed a filter, and a sample with a missing genotype (`./.`) at a site that itself failed a filter.
## How failed calls are counted

What this means for a sample in `fail_bitmap`:

- Its alt alleles are **not** counted in AC, so it does not raise AF.
- It is still an eligible sample, so it **is** counted in AN. The failed call lowers AF by sitting in the denominator without contributing to the numerator.
- Its count is exposed separately as `N_FAIL`.

---
A failed call is judged per **position**, not per allele: a sample with a failed call for any ALT allele at a position has no reliable genotype there. For every allele at that position, the sample:

## Database Creation
- adds nothing to **AC**, even if another record at the position called it with `FILTER=PASS`;
- is left out of **AN**, which counts only the alleles of samples with a usable call;
- is counted in **N_FAIL**, which therefore has the same value on every line of a multi-allelic site;
- stays in **n_eligible**, so `N_HET + N_HOM_ALT + N_HOM_REF + N_FAIL + N_NO_COVERAGE` still accounts for every eligible sample (see [Multi-allelic sites](../getting-started/understanding-output.md#multi-allelic-sites) for the one exception).

The `fail_bitmap` is always written, regardless of build options:

```bash
afquery create-db --manifest manifest.tsv --output-dir ./db/ --genome-build GRCh38
```

---

## Querying N_FAIL

### CLI

The `N_FAIL` count is shown automatically in query output:
Leaving failed calls out of AN keeps AF an estimate over the genotypes that were actually observed. At every position, AN equals the reference alleles plus the AC of all its ALT alleles.

```bash
afquery query --db ./db/ --locus chr1:925952
```

```
chr1:925952 G>A AC=142 AN=2742 AF=0.0518 n_eligible=1371 N_HET=138 N_HOM_ALT=2 N_HOM_REF=1224 N_FAIL=7
chr1:925952 G>A AC=142 AN=2728 AF=0.0521 n_eligible=1371 N_HET=138 N_HOM_ALT=2 N_HOM_REF=1224 N_FAIL=7 N_NO_COVERAGE=0
```

`N_FAIL=7` means 7 eligible samples had a call with FILTER≠PASS at this site. They are part of `n_eligible` (and of AN), so the genotype counts still add up: 138 + 2 + 1224 + 7 = 1371.
The 7 samples with a failed call are part of `n_eligible` (138 + 2 + 1224 + 7 = 1371) but not of AN: 2 × (1371 − 7) = 2728.

### Python API
### When every eligible sample failed

```python
results = db.query("chr1", pos=925952)
for r in results:
print(f"AC={r.AC} AN={r.AN} FAIL={r.N_FAIL}")
if r.N_FAIL > 0:
print(f" Warning: {r.N_FAIL} samples have low-quality calls at this site")
```
If all eligible samples have a failed call at a position, AN is 0 and AF is undefined: `NA` in text and TSV output, `null` in JSON and the Python API, an empty cell in `dump`, and `.` in `AFQUERY_AF`. `N_FAIL` still reports how many samples failed.

`N_FAIL` is always an `int` (default `0`).
---

### Identifying specific FAIL samples
## Counting failed calls as reference

To see which individual samples have FAIL status at a position, use `variant-info`:
`--fail-as-ref` (`fail_as_ref=True` in the Python API) keeps samples with a failed call in AN, counting them as homozygous reference. It is available on `query`, `annotate` and `dump`, and changes only AN and AF: the failed samples still add nothing to AC and are still reported in `N_FAIL`.

```bash
afquery variant-info --db ./db/ --locus chr1:925952
afquery query --db ./db/ --locus chr1:925952 --fail-as-ref
```

Each carrier row shows its `filter` column as `PASS` or `FAIL`, along with sample metadata (technology, phenotype codes). This helps pinpoint whether failures cluster in a specific technology or sample group. See [Variant Info](../guides/variant-info.md) for full options.
```
chr1:925952 G>A AC=142 AN=2742 AF=0.0518 n_eligible=1371 N_HET=138 N_HOM_ALT=2 N_HOM_REF=1224 N_FAIL=7 N_NO_COVERAGE=0
```

---
Because every failed sample is then treated as a non-carrier, the AF is never higher than the default one, which makes it a conservative lower estimate when a site has many failures. AN then no longer equals the reference alleles plus the AC of all ALT alleles.

## VCF Annotation
---

Among the INFO fields `afquery annotate` writes, one reports the failed-call count:
## Finding the failed samples

| Field | Type | Description |
|-------|------|-------------|
| `AFQUERY_N_FAIL` | Integer | Eligible samples with `FILTER≠PASS` at this variant |
`afquery variant-info` lists each sample with a failed call at the position as genotype `alt` with filter `FAIL`, along with its technology and phenotypes. This shows whether failures cluster in one technology or sample group. See [Variant Info](../guides/variant-info.md).

```bash
afquery annotate --db ./db/ --input variants.vcf --output annotated.vcf
afquery variant-info --db ./db/ --locus chr1:925952
```

See [Understanding Output](../getting-started/understanding-output.md#vcf-annotation-fields) for the full set of `AFQUERY_*` fields.

---

## PASS-Only Enforcement
## VCF annotation

PASS-only counting is always enforced and cannot be turned off. It applies to the numerator only: a failed call never adds an alt allele to AC, but the sample stays in the eligible set and so remains in AN. The result is a deliberately conservative AF — a failed carrier weighs on the denominator without lifting the numerator, which is the safe direction for clinical and research use.
`afquery annotate` writes the same per-position values: `AFQUERY_AN` excludes failed calls and `AFQUERY_N_FAIL` counts each failed sample once per site. See [Understanding Output](../getting-started/understanding-output.md#vcf-annotation-fields) for the full set of `AFQUERY_*` fields.

---

Expand Down
15 changes: 12 additions & 3 deletions docs/advanced/ploidy-and-sex-chroms.md
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@ afquery query --db ./db/ --locus chrY:2787758 --sex female
# chrY:2787758 — no results (AN=0 for all variants)

afquery query --db ./db/ --locus chrY:2787758 --sex male
# chrY:2787758 C>T AC=3 AN=856 AF=0.0035 n_eligible=856 N_HET=0 N_HOM_ALT=3 N_HOM_REF=853 N_FAIL=0
# chrY:2787758 C>T AC=3 AN=856 AF=0.0035 n_eligible=856 N_HET=0 N_HOM_ALT=3 N_HOM_REF=853 N_FAIL=0 N_NO_COVERAGE=0
```

### chrX non-PAR
Expand Down Expand Up @@ -107,10 +107,10 @@ 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. 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).
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 have no failed call there). 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.
N_HET, N_HOM_ALT, N_HOM_REF, N_FAIL, and N_NO_COVERAGE are mutually exclusive. A sample with a failed call (FILTER≠PASS) for any allele at the position is counted in N_FAIL only — it does not appear in N_HET, N_HOM_ALT or N_HOM_REF, and is left out of AC and AN.

### chrX non-PAR

Expand Down Expand Up @@ -140,6 +140,15 @@ chrM is haploid for all samples (both sexes contribute AN=1):
!!! warning "N_HET is always 0 on haploid regions"
On chrY, chrM, and chrX non-PAR for males, all carriers are counted in N_HOM_ALT. N_HET is 0 because haploid samples have only one allele copy — there is no heterozygous state. This is correct behavior, not a bug.

### Genotypes inconsistent with ploidy

Ploidy comes from each sample's sex and the position, not from how the VCF writes the genotype: a male carrier on chrX non-PAR counts as one allele whether its GT is `1`, `0/1` or `1/1`. Calls that carry more alleles than the sample has cannot be counted that way:

- on a haploid position, two different ALT alleles (e.g. GT=`1/2` on chrM, or on chrX non-PAR in a male);
- on a diploid position, more than two copies across the records at that position (e.g. `1/1` for one ALT and `0/1` for another).

`create-db` and `update-db --add-samples` record every call of the sample at such a position as a failed call, so it counts in `N_FAIL` and is left out of AC and AN, and log a warning naming the sample and the positions. Running `bcftools +fixploidy` during [preprocessing](../getting-started/preprocessing.md) avoids most of these cases on chrX.

---

## Sex Filter Interaction
Expand Down
2 changes: 1 addition & 1 deletion docs/faq.md
Original file line number Diff line number Diff line change
Expand Up @@ -153,7 +153,7 @@ On chrX non-PAR positions, males contribute AN=1 and females contribute AN=2. Wh

## What does N_FAIL mean in query output?

`N_FAIL` is the count of eligible samples whose call had `FILTER≠PASS` at the position. Shown as `N_FAIL=N` in text output. These samples are excluded from AC, but they stay eligible and still count toward AN. See [FILTER=PASS Tracking](advanced/filter-pass-tracking.md).
`N_FAIL` is the count of eligible samples whose call had `FILTER≠PASS` at the position. Shown as `N_FAIL=N` in text output. These samples are left out of AC and AN, because their genotype there is unknown. See [FILTER=PASS Tracking](advanced/filter-pass-tracking.md).

---

Expand Down
10 changes: 5 additions & 5 deletions docs/getting-started/concepts.md
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ For each variant, AFQuery stores a compressed bitset recording which samples car
```
Samples → Metadata → Filters → Eligible Bitmap
Variants → Genotypes → het/hom Bitmaps
Query: (eligible bitmap) AND (carrier bitmap) → AC; count(eligible) → AN; AC/AN → AF
Query: (eligible bitmap) AND (carrier bitmap) → AC; count(eligible − failed) → AN; AC/AN → AF
```

```mermaid
Expand All @@ -25,7 +25,7 @@ flowchart LR
end

subgraph Result
R["AC = popcount(eligible AND carrier)<br/>AN = popcount(eligible) × ploidy<br/>AF = AC / AN"]
R["AC = popcount(eligible AND carrier)<br/>AN = popcount(eligible − failed) × ploidy<br/>AF = AC / AN"]
end

S --> E
Expand Down Expand Up @@ -55,7 +55,7 @@ For a detailed discussion of the methodological gaps AFQuery addresses, includin
For a given variant (chromosome, position, ref, alt):

- **AC** (Allele Count) — number of copies of the alt allele observed across eligible samples
- **AN** (Allele Number) — total number of alleles examined (2× diploid samples, 1× haploid)
- **AN** (Allele Number) — total number of alleles examined (2× diploid samples, 1× haploid). Samples with a failed call (`FILTER≠PASS`) at the position are left out, as their genotype there is unknown; they are reported as `N_FAIL` (see [FILTER=PASS Tracking](../advanced/filter-pass-tracking.md))
- **AF** (Allele Frequency) — `AC / AN` (None when AN = 0)

A heterozygous carrier contributes AC=1; a homozygous alt carrier contributes AC=2.
Expand Down Expand Up @@ -130,7 +130,7 @@ For each variant row, AFQuery stores per-sample genotype information as [Roaring

- **`het_bitmap`** — sample is heterozygous (`GT=0/1`) with `FILTER=PASS`
- **`hom_bitmap`** — sample is homozygous alt (`GT=1/1`, or `GT=1` on haploid regions) with `FILTER=PASS`
- **`fail_bitmap`** — the sample's call at this site has `FILTER≠PASS`
- **`fail_bitmap`** — the sample's call at this site failed (`FILTER≠PASS`, or a genotype inconsistent with its ploidy); see [FILTER=PASS Tracking](../advanced/filter-pass-tracking.md)

Each sample has a stable integer ID (0-indexed). The bit position in the bitmap equals the sample ID. Databases built with coverage-quality filters carry two additional bitmaps — see [Data Model](../reference/data-model.md#parquet-schema).

Expand Down Expand Up @@ -214,7 +214,7 @@ When planning phenotype codes, consider:

---

PASS-only ingestion is always enforced. See [FILTER=PASS Tracking](../advanced/filter-pass-tracking.md) for details.
Only `FILTER=PASS` calls count towards AC. Failed calls are stored too, and their samples are reported as `N_FAIL` and left out of AN. See [FILTER=PASS Tracking](../advanced/filter-pass-tracking.md) for details.

---

Expand Down
Loading
Loading