From eba2e535d335049c63ae4633294b0384886d1a1c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20L=C3=B3pez=20L=C3=B3pez?= Date: Wed, 14 Oct 2026 12:31:51 +0200 Subject: [PATCH 01/10] fix(query): judge --min-covered evidence per position at multi-allelic sites The build stores filtered_bitmap per ALT allele, counting only that allele's quality-passing carriers. At a multi-allelic site a panel with carriers of one allele but none of another had its non-carriers reported as N_NO_COVERAGE for the second allele and as N_HOM_REF for the first. At positions holding several alleles the --min-covered gate is now recomputed at query time from the pooled quality_pass bitmaps, so every allele of the site gets the same verdict. Single-allele positions keep using the stored bitmap, which is equivalent there. --- docs/advanced/coverage-evidence.md | 2 +- src/afquery/query.py | 14 ++++++++++--- tests/test_multiallelic.py | 32 ++++++++++++++++++++++++++++++ tests/test_variant_info.py | 9 ++++++--- 4 files changed, 50 insertions(+), 7 deletions(-) diff --git a/docs/advanced/coverage-evidence.md b/docs/advanced/coverage-evidence.md index 9e94460..3e8a620 100644 --- a/docs/advanced/coverage-evidence.md +++ b/docs/advanced/coverage-evidence.md @@ -102,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`. 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 diff --git a/src/afquery/query.py b/src/afquery/query.py index c8ae6df..1e7c311 100644 --- a/src/afquery/query.py +++ b/src/afquery/query.py @@ -257,10 +257,18 @@ def _compute_no_coverage_bm( if pass_count < min_pass or observed_count < min_observed: no_cov |= tech_eligible - (all_carriers & tech_eligible) - # Phase 2: stored filtered_bitmap (built per allele, so it may hold - # carriers of another allele at the same position) + # Phase 2: build-time --min-covered. The stored filtered_bitmap counts the + # quality-passing carriers of its own allele only, so at a multi-allelic + # site the gate is recomputed from the carriers of every allele there. if filtered_bm is not None: - no_cov |= (filtered_bm & eligible) - all_carriers + min_covered = self._coverage_filter.get("min_covered", 0) + if site is None or site.quality_pass_bm is None: + no_cov |= (filtered_bm & eligible) - all_carriers + elif min_covered > 0: + for tech_id in self._targeted_capture: + tech_bm = self._tech_bitmaps.get(str(tech_id), BitMap()) + if len(site.quality_pass_bm & tech_bm) < min_covered: + no_cov |= (tech_bm & eligible) - all_carriers # Phase 2: quality_pass gate (--min-quality-evidence K) if quality_pass_bm is not None and min_quality_evidence > 0: diff --git a/tests/test_multiallelic.py b/tests/test_multiallelic.py index 79344b6..abc3b64 100644 --- a/tests/test_multiallelic.py +++ b/tests/test_multiallelic.py @@ -151,3 +151,35 @@ def test_min_pass_failure_never_marks_carriers_uncovered(panel_db): carriers = panel_db.variant_info(*SITE, alt="A", min_pass=2) no_cov = sorted(c.sample_name for c in carriers if c.genotype == "no_coverage") assert no_cov == ["P_R1", "P_R2"] + + +def test_min_covered_judged_per_position(tmp_path): + """--min-covered counts quality-passing carriers of every allele at the + position: the panel carries A, so it was sequenced here for T as well.""" + db = _build_db(tmp_path, [ + ("W_T", "female", "WGS", [("chr1", 5000, "G", "T", "0/1")]), + ("W_R", "female", "WGS", [("chr1", 100, "C", "G", "0/1")]), + ("P_A1", "female", "PANEL", [("chr1", 5000, "G", "A", "0/1")]), + ("P_A2", "male", "PANEL", [("chr1", 5000, "G", "A", "0/1")]), + ("P_R1", "female", "PANEL", [("chr1", 100, "C", "G", "0/1")]), + ("P_R2", "male", "PANEL", [("chr1", 100, "C", "G", "0/1")]), + ], beds={"PANEL": "chr1\t0\t6000\n"}, min_covered=1) + by_alt = {r.variant.alt: r for r in db.query(*SITE)} + for alt in ("A", "T"): + assert by_alt[alt].N_NO_COVERAGE == 0, alt + assert by_alt["A"].N_HOM_REF == 3 # W_R, P_R1, P_R2 + assert by_alt["T"].N_HOM_REF == 3 + + +def test_min_covered_still_flags_unevidenced_panel(tmp_path): + """With no quality-passing carrier of any allele in the panel, its + non-carriers stay uncovered for every allele at the position.""" + db = _build_db(tmp_path, [ + ("W_A", "female", "WGS", [("chr1", 5000, "G", "A", "0/1")]), + ("W_T", "female", "WGS", [("chr1", 5000, "G", "T", "0/1")]), + ("P_R1", "female", "PANEL", [("chr1", 100, "C", "G", "0/1")]), + ("P_R2", "male", "PANEL", [("chr1", 100, "C", "G", "0/1")]), + ], beds={"PANEL": "chr1\t0\t6000\n"}, min_covered=1) + by_alt = {r.variant.alt: r for r in db.query(*SITE)} + for alt in ("A", "T"): + assert by_alt[alt].N_NO_COVERAGE == 2, alt diff --git a/tests/test_variant_info.py b/tests/test_variant_info.py index fc6bea9..3c71661 100644 --- a/tests/test_variant_info.py +++ b/tests/test_variant_info.py @@ -469,8 +469,11 @@ def _write_vcf(path, sample_name, records): f.write(f"{chrom}\t{pos}\t.\t{ref}\t{alt}\t.\t{flt}\t.\tGT\t{gt}\n") -def _build_db(tmp_path, sample_defs, beds=None): - """sample_defs = [(name, sex, tech, records)]; beds = {tech: bed_text}.""" +def _build_db(tmp_path, sample_defs, beds=None, **build_options): + """sample_defs = [(name, sex, tech, records)]; beds = {tech: bed_text}. + + build_options are passed to run_preprocess (e.g. min_covered=1). + """ from afquery.preprocess import run_preprocess bed_dir = tmp_path / "beds" @@ -488,7 +491,7 @@ def _build_db(tmp_path, sample_defs, beds=None): db_path.mkdir() run_preprocess( manifest_path=str(manifest), output_dir=str(db_path), - genome_build="GRCh37", threads=1, bed_dir=str(bed_dir), + genome_build="GRCh37", threads=1, bed_dir=str(bed_dir), **build_options, ) return Database(str(db_path)) From 5ff45ad93d49081e9eb2210939341dee0357b62e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20L=C3=B3pez=20L=C3=B3pez?= Date: Wed, 14 Oct 2026 12:32:32 +0200 Subject: [PATCH 02/10] fix(annotate): apply coverage evidence to alleles absent from the database An input ALT allele with no stored row got N_HOM_REF = eligible minus carriers of the stored alleles, with N_NO_COVERAGE always 0. Under --min-pass, --min-observed or --min-quality-evidence the same samples were uncovered for the stored alleles of the record and hom-ref for the absent one. When other alleles are stored at the position, the absent allele is now computed like any other, with no carriers and the pooled site evidence. Positions with nothing stored keep every eligible sample as hom-ref. --- src/afquery/annotate.py | 27 ++++++++++++++++++++------- tests/test_multiallelic.py | 17 +++++++++++++++++ 2 files changed, 37 insertions(+), 7 deletions(-) diff --git a/src/afquery/annotate.py b/src/afquery/annotate.py index ccede51..3b479de 100644 --- a/src/afquery/annotate.py +++ b/src/afquery/annotate.py @@ -2,6 +2,7 @@ import os import warnings import duckdb +from pyroaring import BitMap from . import storage from .bitmaps import deserialize @@ -89,6 +90,10 @@ def _compute_chunk_annotations( for pos, group in stored_by_pos.items() } + empty = BitMap() + stored_empty = empty if engine._has_coverage_data else None + no_carriers = (empty, empty, empty, stored_empty, stored_empty) + result: dict[tuple[int, str, str], tuple[int, int, bool, int, int, int, int, int]] = {} for pos, ref, alts in records: eligible, AN = pos_data[pos] @@ -105,14 +110,22 @@ def _compute_chunk_annotations( ) result[key] = (s.AC, AN, True, s.N_FAIL, s.N_HET, s.N_HOM_ALT, s.N_HOM_REF, s.N_NO_COVERAGE) + elif pos in stored_by_pos: + # Allele not in Parquet, but others are stored at this position: + # it has no carriers, and the rest of its tallies come from the + # site, as for the stored alleles. + s = engine._variant_stats( + chrom, pos, eligible, no_carriers, + engine._site_evidence(stored_by_pos[pos]), + sf.min_pass, sf.min_observed, sf.min_quality_evidence, + ) + result[key] = (s.AC, AN, False, s.N_FAIL, s.N_HET, s.N_HOM_ALT, + s.N_HOM_REF, s.N_NO_COVERAGE) else: - # Allele not in Parquet: eligible samples are hom-ref unless they - # carry another allele stored at this position. Phase 1/2 filters - # do not apply because this allele has no carriers to evaluate. - stored = stored_by_pos.get(pos) - n_other = (len(engine._site_evidence(stored).carrier_bm & eligible) - if stored else 0) - result[key] = (0, AN, False, 0, 0, 0, len(eligible) - n_other, 0) + # Nothing stored at this position: every eligible sample is + # hom-ref. Coverage-evidence filters do not apply, as there is + # no call to evaluate. + result[key] = (0, AN, False, 0, 0, 0, len(eligible), 0) return result diff --git a/tests/test_multiallelic.py b/tests/test_multiallelic.py index abc3b64..369b021 100644 --- a/tests/test_multiallelic.py +++ b/tests/test_multiallelic.py @@ -183,3 +183,20 @@ def test_min_covered_still_flags_unevidenced_panel(tmp_path): by_alt = {r.variant.alt: r for r in db.query(*SITE)} for alt in ("A", "T"): assert by_alt[alt].N_NO_COVERAGE == 2, alt + + +def test_annotate_absent_allele_shares_site_coverage(panel_db, tmp_path): + """C is not stored: under --min-pass 2 the panel's non-carriers are as + uncovered for C as they are for the stored A and T.""" + vcf_in = tmp_path / "in.vcf" + vcf_in.write_text( + "##fileformat=VCFv4.2\n" + "##contig=\n" + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n" + "chr1\t5000\t.\tG\tA,T,C\t.\tPASS\t.\n" + ) + out = tmp_path / "out.vcf" + panel_db.annotate_vcf(str(vcf_in), str(out), min_pass=2) + [v] = list(cyvcf2.VCF(str(out))) + assert tuple(v.INFO.get("AFQUERY_N_NO_COVERAGE")) == (2, 2, 2) + assert tuple(v.INFO.get("AFQUERY_N_HOM_REF")) == (2, 2, 2) From 2d74b781f4e8b5a36b4b805ab777442ff2c57944 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20L=C3=B3pez=20L=C3=B3pez?= Date: Wed, 14 Oct 2026 12:47:56 +0200 Subject: [PATCH 03/10] feat(query)!: exclude failed calls from AN and count N_FAIL per site A sample with a failed call (FILTER!=PASS) at a position was left out of AC but counted as reference in AN, so AF was computed as if every failed sample were hom-ref. At multi-allelic sites the alleles did not add up either: AN exceeded the reference alleles plus the AC of all ALT alleles by the alleles of the failed samples. A failed call for any allele at a position now makes the sample's genotype there unknown, for every allele of the site: it is counted in N_FAIL, which is the same on every line of the position, and left out of AC and AN. AN therefore equals the reference alleles plus the AC of all ALT alleles. When every eligible sample failed, AN is 0 and AF is undefined (NA in query output, an empty cell in dump, "." in annotate). annotate wrote AFQUERY_N_FAIL as the sum of the per-allele counts, so a sample masked as ./. at a failed multi-allelic site was counted once per allele. It now reports the per-site N_FAIL. variant-info lists the failed samples for every allele of the position, matching N_FAIL. The test fixture stored some samples in both het and fail bitmaps, which a real build never does; it and its source VCFs now agree. BREAKING CHANGE: AN no longer includes samples with a failed call, so AN and AF change at every position where N_FAIL > 0. N_FAIL counts failed samples per position rather than per allele. --- README.md | 2 +- docs/advanced/coverage-evidence.md | 5 +- docs/advanced/debugging-results.md | 4 +- docs/advanced/filter-pass-tracking.md | 90 +++++--------------- docs/advanced/ploidy-and-sex-chroms.md | 4 +- docs/faq.md | 2 +- docs/getting-started/concepts.md | 8 +- docs/getting-started/preprocessing.md | 2 +- docs/getting-started/quickstart.md | 2 +- docs/getting-started/tutorial.md | 37 ++++---- docs/getting-started/understanding-output.md | 26 ++++-- docs/guides/annotate-vcf.md | 8 +- docs/guides/create-database.md | 2 +- docs/guides/dump-export.md | 2 + docs/guides/query.md | 2 +- docs/guides/variant-info.md | 2 +- docs/reference/data-model.md | 2 +- docs/reference/glossary.md | 8 +- docs/reference/python-api.md | 4 +- docs/use-cases/acmg-use-cases.md | 2 +- src/afquery/annotate.py | 38 +++++---- src/afquery/dump.py | 21 +++-- src/afquery/query.py | 69 +++++++++------ tests/conftest.py | 9 +- tests/data/expected_results.json | 8 +- tests/data/vcfs/S00.vcf | 3 +- tests/data/vcfs/S01.vcf | 3 +- tests/data/vcfs/S03.vcf | 3 +- tests/data/vcfs/S05.vcf | 3 +- tests/oracle.py | 26 ++++-- tests/test_annotate.py | 25 +++--- tests/test_cli.py | 16 ++-- tests/test_dump.py | 25 ++++-- tests/test_haploid_stats.py | 12 +-- tests/test_multiallelic.py | 83 ++++++++++++++++++ tests/test_no_coverage.py | 3 +- tests/test_oracle_consistency.py | 20 +++-- tests/test_sample_filter.py | 12 +-- tests/test_update.py | 13 ++- tests/test_variant_info.py | 21 +++-- 40 files changed, 368 insertions(+), 259 deletions(-) diff --git a/README.md b/README.md index 7bee2a0..8c9aa00 100644 --- a/README.md +++ b/README.md @@ -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`). diff --git a/docs/advanced/coverage-evidence.md b/docs/advanced/coverage-evidence.md index 3e8a620..c4cb7ec 100644 --- a/docs/advanced/coverage-evidence.md +++ b/docs/advanced/coverage-evidence.md @@ -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 diff --git a/docs/advanced/debugging-results.md b/docs/advanced/debugging-results.md index 03834c4..49e7a13 100644 --- a/docs/advanced/debugging-results.md +++ b/docs/advanced/debugging-results.md @@ -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: diff --git a/docs/advanced/filter-pass-tracking.md b/docs/advanced/filter-pass-tracking.md index 11d03ca..6fff8e6 100644 --- a/docs/advanced/filter-pass-tracking.md +++ b/docs/advanced/filter-pass-tracking.md @@ -1,106 +1,60 @@ # 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 `.`; or +- a missing genotype (`./.`) at a site that failed a filter. -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. - -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`. - ---- - -## Database Creation - -The `fail_bitmap` is always written, regardless of build options: - -```bash -afquery create-db --manifest manifest.tsv --output-dir ./db/ --genome-build GRCh38 -``` - ---- +## How failed calls are counted -## Querying 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: -### CLI +- 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 `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_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. - -### Python API - -```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") -``` +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. -`N_FAIL` is always an `int` (default `0`). +### When every eligible sample failed -### Identifying specific FAIL samples - -To see which individual samples have FAIL status at a position, use `variant-info`: - -```bash -afquery variant-info --db ./db/ --locus chr1:925952 -``` - -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. +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. --- -## VCF Annotation +## Finding the failed samples -Among the INFO fields `afquery annotate` writes, one reports the failed-call count: - -| 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. --- diff --git a/docs/advanced/ploidy-and-sex-chroms.md b/docs/advanced/ploidy-and-sex-chroms.md index 36acbcf..efbbeac 100644 --- a/docs/advanced/ploidy-and-sex-chroms.md +++ b/docs/advanced/ploidy-and-sex-chroms.md @@ -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 diff --git a/docs/faq.md b/docs/faq.md index 0fdab9b..f15e0b3 100644 --- a/docs/faq.md +++ b/docs/faq.md @@ -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). --- diff --git a/docs/getting-started/concepts.md b/docs/getting-started/concepts.md index fb1d5be..0474070 100644 --- a/docs/getting-started/concepts.md +++ b/docs/getting-started/concepts.md @@ -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 @@ -25,7 +25,7 @@ flowchart LR end subgraph Result - R["AC = popcount(eligible AND carrier)
AN = popcount(eligible) × ploidy
AF = AC / AN"] + R["AC = popcount(eligible AND carrier)
AN = popcount(eligible − failed) × ploidy
AF = AC / AN"] end S --> E @@ -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. @@ -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 has `FILTER≠PASS`; it counts in neither AC nor AN 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). diff --git a/docs/getting-started/preprocessing.md b/docs/getting-started/preprocessing.md index 08e6f64..e1f564d 100644 --- a/docs/getting-started/preprocessing.md +++ b/docs/getting-started/preprocessing.md @@ -38,7 +38,7 @@ Male samples should have haploid genotype calls at chrX non-PAR regions (GT=`1`, ### Filtering and stripping -- **Non-PASS genotypes**: Masking them as missing ensures that low-quality calls do not inflate AC. AFQuery tracks these as N_FAIL. +- **Non-PASS genotypes**: Masking every genotype at a failed site as missing ensures that low-quality calls do not inflate AC, and that no sample there is taken for hom-ref. AFQuery tracks these samples as N_FAIL and leaves them out of AN. - **Homozygous reference calls**: Removing ref/ref genotypes reduces file size and speeds ingestion; they contribute AC=0 and are not needed - **INFO fields**: Stripping INFO reduces file size and speeds ingestion. Additionally, malformed or non-standard INFO fields produced by some variant callers can break downstream parsing; stripping them pre-emptively prevents these errors. - **FORMAT fields**: Only `GT`, `DP`, and `GQ` are preserved. `GT` is the genotype required by AFQuery for all queries. `DP` and `GQ` are read by the [coverage-evidence](../advanced/coverage-evidence.md) quality flags (`afquery create-db --min-dp / --min-gq / --min-qual / --min-covered`). VCFs without `DP`/`GQ` are still valid; their carriers simply contribute no quality evidence to the cohort. All other FORMAT fields (`PL`, `AD`, etc.) are dropped to reduce file size. diff --git a/docs/getting-started/quickstart.md b/docs/getting-started/quickstart.md index 1a3ace0..ad416db 100644 --- a/docs/getting-started/quickstart.md +++ b/docs/getting-started/quickstart.md @@ -134,7 +134,7 @@ The output VCF gains INFO fields (see [Annotate a VCF](../guides/annotate-vcf.md | `AFQUERY_N_HET` | A (per ALT) | Heterozygous sample count | | `AFQUERY_N_HOM_ALT` | A (per ALT) | Homozygous alt sample count | | `AFQUERY_N_HOM_REF` | A (per ALT) | Homozygous ref sample count | -| `AFQUERY_N_FAIL` | 1 (per site) | Samples with FILTER≠PASS | +| `AFQUERY_N_FAIL` | 1 (per site) | Samples with a failed call (FILTER≠PASS), left out of AN | | `AFQUERY_N_NO_COVERAGE` | A (per ALT) | Eligible samples lacking coverage evidence (0 unless a coverage-evidence filter is active) | --- diff --git a/docs/getting-started/tutorial.md b/docs/getting-started/tutorial.md index df85551..e3fc64e 100644 --- a/docs/getting-started/tutorial.md +++ b/docs/getting-started/tutorial.md @@ -107,15 +107,10 @@ afquery query --db ./demo_db/ --locus chr1:946000 Example output: ``` -chr1:946000 T>C AC=2 AN=20 AF=0.1000 n_eligible=10 N_HET=2 N_HOM_ALT=0 N_HOM_REF=6 N_FAIL=2 +chr1:946000 T>C AC=2 AN=16 AF=0.1250 n_eligible=10 N_HET=2 N_HOM_ALT=0 N_HOM_REF=6 N_FAIL=2 N_NO_COVERAGE=0 ``` -Here, 2 samples were flagged as `LowQual` in their source VCF at this position. N_FAIL samples are: - -- **Counted in AN** — they are eligible and contribute to the denominator -- **Excluded from AC** — their alleles are not counted in the numerator - -This makes AF a conservative estimate: `AF=2/20=0.1000` even though 2 additional samples may carry the variant. Verify: 10 = 2+0+6+2 ✓. +Here, 2 samples were flagged as `LowQual` in their source VCF at this position. Their genotype there is unknown, so they are left out of both AC and AN: `AN=2×(10−2)=16` and `AF=2/16=0.1250`. They remain eligible, so the accounting identity still holds: 10 = 2+0+6+2+0 ✓. See [FILTER=PASS Tracking](../advanced/filter-pass-tracking.md). !!! tip N_FAIL > 0 is a signal to inspect source VCF quality at this site. A high N_FAIL relative to n_eligible may indicate a systematic sequencing artifact. See [Understanding Output](understanding-output.md) for thresholds. @@ -131,9 +126,9 @@ afquery query --db ./demo_db/ --region chr1:900000-950000 Example output: ``` -chr1:925100 C>T AC=4 AN=20 AF=0.2000 n_eligible=10 N_HET=4 N_HOM_ALT=0 N_HOM_REF=6 N_FAIL=0 -chr1:925952 G>A AC=6 AN=20 AF=0.3000 n_eligible=10 N_HET=4 N_HOM_ALT=1 N_HOM_REF=5 N_FAIL=0 -chr1:946000 T>C AC=2 AN=20 AF=0.1000 n_eligible=10 N_HET=2 N_HOM_ALT=0 N_HOM_REF=6 N_FAIL=2 +chr1:925100 C>T AC=4 AN=20 AF=0.2000 n_eligible=10 N_HET=4 N_HOM_ALT=0 N_HOM_REF=6 N_FAIL=0 N_NO_COVERAGE=0 +chr1:925952 G>A AC=6 AN=20 AF=0.3000 n_eligible=10 N_HET=4 N_HOM_ALT=1 N_HOM_REF=5 N_FAIL=0 N_NO_COVERAGE=0 +chr1:946000 T>C AC=2 AN=16 AF=0.1250 n_eligible=10 N_HET=2 N_HOM_ALT=0 N_HOM_REF=6 N_FAIL=2 N_NO_COVERAGE=0 ``` --- @@ -325,18 +320,18 @@ afquery annotate \ --output ./annotated_demo.vcf ``` -The output VCF gains `AFQUERY_AC`, `AFQUERY_AN`, `AFQUERY_AF`, `AFQUERY_N_HET`, `AFQUERY_N_HOM_ALT`, `AFQUERY_N_HOM_REF`, and `AFQUERY_N_FAIL` INFO fields. Inspect the chr1 variants: +The output VCF gains `AFQUERY_AC`, `AFQUERY_AN`, `AFQUERY_AF`, `AFQUERY_N_HET`, `AFQUERY_N_HOM_ALT`, `AFQUERY_N_HOM_REF`, `AFQUERY_N_FAIL`, and `AFQUERY_N_NO_COVERAGE` INFO fields. Inspect the chr1 variants: ```bash -grep -v "^##" ./annotated_demo.vcf | grep "^chr1" | head -5 +grep -v "^##" ./annotated_demo.vcf | grep "^chr1" | head -3 ``` Example output (columns: CHROM, POS, ID, REF, ALT, QUAL, FILTER, INFO, FORMAT, DEMO_001): ``` -chr1 925100 . C T 50 PASS AFQUERY_AC=4;AFQUERY_AN=20;AFQUERY_AF=0.2;AFQUERY_N_HET=4;AFQUERY_N_HOM_ALT=0;AFQUERY_N_HOM_REF=6;AFQUERY_N_FAIL=0 GT 0/1 -chr1 925952 . G A 50 PASS AFQUERY_AC=6;AFQUERY_AN=20;AFQUERY_AF=0.3;AFQUERY_N_HET=4;AFQUERY_N_HOM_ALT=1;AFQUERY_N_HOM_REF=5;AFQUERY_N_FAIL=0 GT 0/1 -chr1 946000 . T C 50 PASS AFQUERY_AC=2;AFQUERY_AN=20;AFQUERY_AF=0.1;AFQUERY_N_HET=2;AFQUERY_N_HOM_ALT=0;AFQUERY_N_HOM_REF=6;AFQUERY_N_FAIL=2 GT 0/1 +chr1 925100 . C T 50 PASS AFQUERY_AC=4;AFQUERY_AN=20;AFQUERY_AF=0.2;AFQUERY_N_HET=4;AFQUERY_N_HOM_ALT=0;AFQUERY_N_HOM_REF=6;AFQUERY_N_FAIL=0;AFQUERY_N_NO_COVERAGE=0 GT 0/1 +chr1 925952 . G A 50 PASS AFQUERY_AC=6;AFQUERY_AN=20;AFQUERY_AF=0.3;AFQUERY_N_HET=4;AFQUERY_N_HOM_ALT=1;AFQUERY_N_HOM_REF=5;AFQUERY_N_FAIL=0;AFQUERY_N_NO_COVERAGE=0 GT 0/1 +chr1 946000 . T C 50 PASS AFQUERY_AC=2;AFQUERY_AN=16;AFQUERY_AF=0.125;AFQUERY_N_HET=2;AFQUERY_N_HOM_ALT=0;AFQUERY_N_HOM_REF=6;AFQUERY_N_FAIL=2;AFQUERY_N_NO_COVERAGE=0 GT 0/1 ``` `AFQUERY_N_FAIL=2` for `chr1:946000` matches the query output from Step 4. The annotated values reflect cohort-wide frequencies across all 10 samples, not just DEMO_001. @@ -362,12 +357,12 @@ head -6 demo_dump.csv Example output: ``` -chrom,pos,ref,alt,AC,AN,AF,N_HET,N_HOM_ALT,N_HOM_REF,N_FAIL -chr1,925100,C,T,4,20,0.2,4,0,6,0 -chr1,925952,G,A,6,20,0.3,4,1,5,0 -chr1,946000,T,C,2,20,0.1,2,0,6,2 -chr1,1225000,A,G,5,20,0.25,5,0,5,0 -chr1,1399914,G,T,2,8,0.25,2,0,2,0 +chrom,pos,ref,alt,AC,AN,AF,N_HET,N_HOM_ALT,N_HOM_REF,N_FAIL,N_NO_COVERAGE +chr1,925100,C,T,4,20,0.2,4,0,6,0,0 +chr1,925952,G,A,6,20,0.3,4,1,5,0,0 +chr1,946000,T,C,2,16,0.125,2,0,6,2,0 +chr1,1225000,A,G,5,20,0.25,5,0,5,0,0 +chr1,1399914,G,T,2,8,0.25,2,0,2,0,0 ``` Note that `chr1:1399914` has AN=8 even in the full-cohort dump — technology-aware AN is applied automatically (WES samples are excluded for positions outside their capture regions). diff --git a/docs/getting-started/understanding-output.md b/docs/getting-started/understanding-output.md index b7d7781..bc74475 100644 --- a/docs/getting-started/understanding-output.md +++ b/docs/getting-started/understanding-output.md @@ -9,14 +9,14 @@ This page explains what each field in AFQuery output means and how to interpret | Field | Type | Description | |-------|------|-------------| | **AC** | int | Allele count — number of alt allele copies in eligible samples | -| **AN** | int | Allele number — total alleles examined (adjusted for ploidy and eligible samples) | -| **AF** | float | Allele frequency — `AC / AN`. `None` when AN=0 | +| **AN** | int | Allele number — alleles of the eligible samples with a usable call at this position, adjusted for ploidy. Samples in `N_FAIL` are left out | +| **AF** | float | Allele frequency — `AC / AN`. `NA` (text/TSV) or `null` (JSON) 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). 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). | +| **N_FAIL** | int | Number of eligible samples with a failed call (FILTER≠PASS) for any allele at this position — the same value for every allele there. These samples are counted *only* in N_FAIL, not in N_HET, N_HOM_ALT or N_HOM_REF, and are left out of AC and AN. See [FILTER=PASS Tracking](../advanced/filter-pass-tracking.md). | +| **N_NO_COVERAGE** | int | Number of eligible samples whose tech lacks coverage evidence at this position. Excluded from `N_HOM_REF`; they still count in AN. Always `0` unless a coverage-evidence filter is active. See [Coverage Evidence](../advanced/coverage-evidence.md). | --- @@ -74,11 +74,14 @@ afquery query --db ./db/ --locus chr1:925952 --ref G --alt A --format json ### AN=0 and AF=None -AN=0 means no eligible samples have coverage at this position. This happens when: +AN=0 means no eligible sample has a usable call at this position. This happens when: - All eligible samples are WES or panels and the position is outside capture regions - The phenotype/sex/technology filter excludes all samples - The chromosome is not in the database +- Every eligible sample has a failed call there (`N_FAIL` = `n_eligible`) + +In the first three cases `query` returns no line for the position; in the last one it reports the variant with AN=0. !!! warning "AN=0 does not mean the variant is absent" @@ -94,6 +97,11 @@ reference for this allele, so it appears in none of `N_HET`, `N_HOM_ALT`, `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. +AN and `N_FAIL` describe the position, so they are the same on every line of a +multi-allelic site, and AN equals the reference alleles plus the AC of all its +ALT alleles. A sample with a failed call for any of the alleles is counted in +`N_FAIL` on every line. + A sample with two different ALT alleles (`GT=1/2`) is counted as heterozygous on both lines. @@ -124,7 +132,7 @@ AC=0 with a high AN (e.g., AN=4000) means the variant was **genuinely not observ ### N_FAIL > 0 -When N_FAIL > 0, some eligible samples had a non-PASS filter at this position in their source VCF. A high N_FAIL relative to AN may indicate a problematic site (e.g., systematic sequencing artifacts). Consider filtering positions with more than 10% of failing (`N_FAIL / (AN / 2) > 0.1`). +When N_FAIL > 0, some eligible samples had a non-PASS filter at this position in their source VCF, and AF is computed over the remaining samples. A high N_FAIL relative to `n_eligible` may indicate a problematic site (e.g., systematic sequencing artifacts). Consider filtering positions where more than 10% of eligible samples failed (`N_FAIL / n_eligible > 0.1`). --- @@ -135,12 +143,12 @@ When using `afquery annotate`, the following INFO fields are added to each varia | INFO field | Number | Description | |-----------|--------|-------------| | `AFQUERY_AC` | A (per ALT) | Allele count — one value per ALT allele | -| `AFQUERY_AN` | 1 (per site) | Allele number — shared across all ALT alleles | -| `AFQUERY_AF` | A (per ALT) | Allele frequency — one value per ALT allele | +| `AFQUERY_AN` | 1 (per site) | Allele number — shared across all ALT alleles; excludes samples in `AFQUERY_N_FAIL` | +| `AFQUERY_AF` | A (per ALT) | Allele frequency — one value per ALT allele; `.` when `AFQUERY_AN=0` | | `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, 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_FAIL` | 1 (per site) | Samples with a failed call for any ALT allele at the position, each counted once | | `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). | !!! note "Multi-allelic sites" diff --git a/docs/guides/annotate-vcf.md b/docs/guides/annotate-vcf.md index e9d3b4a..f3c1f89 100644 --- a/docs/guides/annotate-vcf.md +++ b/docs/guides/annotate-vcf.md @@ -21,13 +21,13 @@ afquery annotate \ | Field | Type | Number | Description | |-------|------|--------|-------------| | `AFQUERY_AC` | Integer | A (per ALT) | Allele count in eligible samples | -| `AFQUERY_AN` | Integer | 1 (per site) | Allele number (total alleles examined) | -| `AFQUERY_AF` | Float | A (per ALT) | Allele frequency (`AC / AN`) | +| `AFQUERY_AN` | Integer | 1 (per site) | Allele number (total alleles examined, excluding samples in `AFQUERY_N_FAIL`) | +| `AFQUERY_AF` | Float | A (per ALT) | Allele frequency (`AC / AN`); `.` when `AFQUERY_AN=0` | | `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, 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). | +| `AFQUERY_N_FAIL` | Integer | 1 (per site) | Eligible samples with a failed call (FILTER≠PASS) for any ALT allele at the position, each counted once. Excluded from AC and 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`; still counted in AN. Always `0` unless a coverage-evidence filter is active. See [Coverage Evidence](../advanced/coverage-evidence.md). | !!! note "Multi-allelic sites" Number=A fields have one value per ALT allele (comma-separated for multi-allelic sites). Number=1 fields are shared across all ALT alleles at the same position. diff --git a/docs/guides/create-database.md b/docs/guides/create-database.md index 6e48a08..a442090 100644 --- a/docs/guides/create-database.md +++ b/docs/guides/create-database.md @@ -93,7 +93,7 @@ afquery create-db --manifest manifest.tsv --output-dir ./db/ --genome-build GRCh ## FILTER=PASS Behavior -Only `FILTER=PASS` calls (or calls with no FILTER field) contribute to AC, and therefore to AF. AN is not affected — it counts every eligible sample, failed calls included. Calls that fail a filter are tracked in `fail_bitmap` and surfaced as `N_FAIL`. PASS-only counting is always enforced — there is currently no CLI option to change this behaviour. +Only `FILTER=PASS` calls (or calls with no FILTER field) contribute to AC. Calls that fail a filter are stored in `fail_bitmap`; at query time those samples are left out of AC and AN and reported as `N_FAIL`. See [FILTER=PASS Tracking](../advanced/filter-pass-tracking.md) for details. diff --git a/docs/guides/dump-export.md b/docs/guides/dump-export.md index c8ddf10..2368609 100644 --- a/docs/guides/dump-export.md +++ b/docs/guides/dump-export.md @@ -70,6 +70,8 @@ chrom,pos,ref,alt,AC,AN,AF,N_HET,N_HOM_ALT,N_HOM_REF,N_FAIL,N_NO_COVERAGE `N_NO_COVERAGE` is always emitted but is `0` unless a coverage-evidence filter is active (see [Coverage Evidence](../advanced/coverage-evidence.md)). +`AF` and every `AF_` column are left empty when the matching `AN` is 0, +i.e. when no eligible sample of that group has a usable call at the position. === "--by-sex" diff --git a/docs/guides/query.md b/docs/guides/query.md index 5b7d064..3f83dc3 100644 --- a/docs/guides/query.md +++ b/docs/guides/query.md @@ -151,7 +151,7 @@ hom-ref. With standard variant-only VCFs that assumption can be wrong: a missing position may simply mean the sample was not sequenced deeply enough at that locus. Three optional flags let you trade hom-ref aggressiveness for confidence. Samples that fall below a threshold are reported in **N_NO_COVERAGE** instead of N_HOM_REF -(they remain in `eligible` and `AN`, like `N_FAIL`). +(they remain in `eligible` and `AN`). | Flag | Meaning | |------|---------| diff --git a/docs/guides/variant-info.md b/docs/guides/variant-info.md index f0953ae..cea5429 100644 --- a/docs/guides/variant-info.md +++ b/docs/guides/variant-info.md @@ -141,7 +141,7 @@ When `--ref` and `--alt` are specified, the `variant` block contains the actual |---|---| | `het` | Heterozygous carrier, FILTER=PASS | | `hom` | Homozygous alt carrier, FILTER=PASS | -| `alt` | Non-ref carrier with FILTER≠PASS (ploidy unknown) | +| `alt` | Failed call (FILTER≠PASS) for any allele at this position; genotype unknown. Listed for every allele of the position, as in `N_FAIL` | | `no_coverage` | Sample's tech lacks coverage evidence at this position; not a carrier. Only appears when a coverage-evidence filter (`--min-pass`, `--min-observed`, `--min-quality-evidence`, or build-time `--min-covered`) is active. The FILTER column is empty (text/tsv) or `null` (JSON) — `PASS`/`FAIL` does not apply because there is no call. See [Coverage Evidence](../advanced/coverage-evidence.md). | --- diff --git a/docs/reference/data-model.md b/docs/reference/data-model.md index d726731..5ea68fd 100644 --- a/docs/reference/data-model.md +++ b/docs/reference/data-model.md @@ -136,7 +136,7 @@ Bitmaps use the [Roaring Bitmap](https://roaringbitmap.org/) format, serialized - Bit position = sample ID (0-indexed integer) - `het_bitmap`: bit set iff the sample is heterozygous with `FILTER=PASS` - `hom_bitmap`: bit set iff the sample is homozygous alt with `FILTER=PASS` -- `fail_bitmap`: bit set iff the sample's call has `FILTER≠PASS` — either a carrier whose call failed a filter, or a missing genotype (`./.`) at a failed site +- `fail_bitmap`: bit set iff the sample's call has `FILTER≠PASS` — either a carrier whose call failed a filter, or a missing genotype (`./.`) at a failed site. A sample set in the `fail_bitmap` of any allele at a position is left out of AC and AN for every allele there - `filtered_bitmap`, `quality_pass_bitmap`: see the Parquet schema above; both are empty unless the database was built with coverage-quality filters To deserialize in Python: diff --git a/docs/reference/glossary.md b/docs/reference/glossary.md index 0db3816..461af3b 100644 --- a/docs/reference/glossary.md +++ b/docs/reference/glossary.md @@ -10,7 +10,7 @@ Ratio of AC to AN (`AC / AN`). Represents the proportion of alleles that are the ## AN (Allele Number) -Total number of alleles examined across eligible samples. For autosomes, AN = 2 × number of eligible diploid samples. For sex chromosomes, AN depends on ploidy. See [Ploidy & Sex Chromosomes](../advanced/ploidy-and-sex-chroms.md). +Total number of alleles examined across eligible samples, leaving out samples with a failed call at the position ([N_FAIL](#n_fail)). For autosomes, AN = 2 × (eligible samples − N_FAIL). For sex chromosomes, AN depends on ploidy. See [Ploidy & Sex Chromosomes](../advanced/ploidy-and-sex-chroms.md). ## Bitmap @@ -30,15 +30,15 @@ An interval tree (stored as a `.pkl` file) built from a BED file that defines th ## Eligible Sample -A sample that passes all query filters (sex, phenotype, technology) and has coverage at the queried position. Only eligible samples contribute to AC and AN. See [Sample Filtering](../guides/sample-filtering.md). +A sample that passes all query filters (sex, phenotype, technology) and has coverage at the queried position. Only eligible samples contribute to AC and AN, and only if their call there did not fail. See [Sample Filtering](../guides/sample-filtering.md). ## N_FAIL -The count of eligible samples whose call at a position had FILTER≠PASS in the source VCF. These samples are excluded from AC, but they remain eligible and so still count toward AN. See [FILTER=PASS Tracking](../advanced/filter-pass-tracking.md). +The count of eligible samples with a failed call (FILTER≠PASS) for any allele at a position. Their genotype there is unknown, so they are left out of AC and AN for every allele at the position. See [FILTER=PASS Tracking](../advanced/filter-pass-tracking.md). ## N_NO_COVERAGE -The count 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 (`--min-pass`, `--min-observed`, `--min-quality-evidence`, or build-time `--min-covered`) is active. See [Coverage Evidence](../advanced/coverage-evidence.md). +The count of eligible samples whose tech lacks coverage evidence at this position. Excluded from `N_HOM_REF`; they still count in AN. Always `0` unless a coverage-evidence filter (`--min-pass`, `--min-observed`, `--min-quality-evidence`, or build-time `--min-covered`) is active. See [Coverage Evidence](../advanced/coverage-evidence.md). ## Manifest diff --git a/docs/reference/python-api.md b/docs/reference/python-api.md index f4e2840..ba31a06 100644 --- a/docs/reference/python-api.md +++ b/docs/reference/python-api.md @@ -447,13 +447,13 @@ Return all distinct phenotype codes present in the database. class QueryResult: variant: VariantKey # chrom, pos, ref, alt AC: int # Allele count - AN: int # Allele number + AN: int # Allele number (samples in N_FAIL left out) AF: float | None # Allele frequency (None if AN == 0) n_samples_eligible: int # Number of eligible samples at this position N_HET: int # Heterozygous count N_HOM_ALT: int # Homozygous alt count N_HOM_REF: int # Homozygous ref count - N_FAIL: int # Eligible samples whose call had FILTER≠PASS (excluded from AC, kept in AN) + N_FAIL: int # Eligible samples with a failed call at the position (out of AC and AN) N_NO_COVERAGE: int # Eligible samples whose tech lacks evidence (excluded from N_HOM_REF) ``` diff --git a/docs/use-cases/acmg-use-cases.md b/docs/use-cases/acmg-use-cases.md index 5ebb3b7..0f183bd 100644 --- a/docs/use-cases/acmg-use-cases.md +++ b/docs/use-cases/acmg-use-cases.md @@ -207,7 +207,7 @@ On chrX non-PAR regions, males contribute AN=1 and females AN=2. A cohort with 8 ### High N_FAIL at the Variant Site -`N_FAIL` counts eligible samples whose call at the site failed quality filters (`FILTER≠PASS`). These samples are excluded from AC but stay in AN, so they cannot inflate AF — if anything they make it slightly more conservative. A high N_FAIL still indicates the site has systematic QC problems, and warrants caution before applying ACMG criteria. +`N_FAIL` counts eligible samples whose call at the site failed quality filters (`FILTER≠PASS`). These samples are left out of both AC and AN, so AF is computed over the samples with a usable call. A high N_FAIL means that AF rests on fewer samples and that the site may have systematic QC problems, which warrants caution before applying ACMG criteria. | N_FAIL / n_eligible | Interpretation | |---|---| diff --git a/src/afquery/annotate.py b/src/afquery/annotate.py index 3b479de..d0359c2 100644 --- a/src/afquery/annotate.py +++ b/src/afquery/annotate.py @@ -96,19 +96,19 @@ def _compute_chunk_annotations( result: dict[tuple[int, str, str], tuple[int, int, bool, int, int, int, int, int]] = {} for pos, ref, alts in records: - eligible, AN = pos_data[pos] + eligible, eligible_AN = pos_data[pos] for alt in alts: key = (pos, ref, alt) if key in result: continue # dedup - if AN == 0: + if eligible_AN == 0: result[key] = (0, 0, False, 0, 0, 0, 0, 0) elif key in unpacked: s = engine._variant_stats( chrom, pos, eligible, unpacked[key], sites[pos], sf.min_pass, sf.min_observed, sf.min_quality_evidence, ) - result[key] = (s.AC, AN, True, s.N_FAIL, s.N_HET, s.N_HOM_ALT, + result[key] = (s.AC, s.AN, True, s.N_FAIL, s.N_HET, s.N_HOM_ALT, s.N_HOM_REF, s.N_NO_COVERAGE) elif pos in stored_by_pos: # Allele not in Parquet, but others are stored at this position: @@ -119,13 +119,13 @@ def _compute_chunk_annotations( engine._site_evidence(stored_by_pos[pos]), sf.min_pass, sf.min_observed, sf.min_quality_evidence, ) - result[key] = (s.AC, AN, False, s.N_FAIL, s.N_HET, s.N_HOM_ALT, + result[key] = (s.AC, s.AN, False, s.N_FAIL, s.N_HET, s.N_HOM_ALT, s.N_HOM_REF, s.N_NO_COVERAGE) else: # Nothing stored at this position: every eligible sample is # hom-ref. Coverage-evidence filters do not apply, as there is # no call to evaluate. - result[key] = (0, AN, False, 0, 0, 0, len(eligible), 0) + result[key] = (0, eligible_AN, False, 0, 0, 0, len(eligible), 0) return result @@ -156,7 +156,8 @@ def annotate_vcf( }) vcf.add_info_to_header({ "ID": "AFQUERY_AN", "Number": "1", "Type": "Integer", - "Description": "Allele number in eligible samples", + "Description": "Allele number in eligible samples, excluding samples with a " + "failed call at the position", }) vcf.add_info_to_header({ "ID": "AFQUERY_AF", "Number": "A", "Type": "Float", @@ -177,7 +178,8 @@ def annotate_vcf( }) vcf.add_info_to_header({ "ID": "AFQUERY_N_FAIL", "Number": "1", "Type": "Integer", - "Description": "Eligible samples with FILTER!=PASS for this variant", + "Description": "Eligible samples with a failed call (FILTER!=PASS) for any " + "allele at the position", }) vcf.add_info_to_header({ "ID": "AFQUERY_N_NO_COVERAGE", "Number": "A", "Type": "Integer", @@ -241,12 +243,14 @@ def annotate_vcf( pos = variant.POS stats["n_variants"] += 1 - # AN is per-position; read from any alt's entry (same value for all alts) - AN = 0 + # AN and N_FAIL are per position: every alt's entry holds the same + # values. AN is 0 with N_FAIL 0 only when no sample is eligible. + AN = N_FAIL = 0 if variant.ALT: - AN = ann.get((pos, variant.REF, variant.ALT[0]), (0, 0, False))[1] + site_entry = ann.get((pos, variant.REF, variant.ALT[0]), (0, 0, False, 0)) + AN, N_FAIL = site_entry[1], site_entry[3] - if AN == 0: + if AN == 0 and N_FAIL == 0: stats["n_uncovered"] += 1 variant.INFO["AFQUERY_AN"] = 0 writer.write_record(variant) @@ -259,31 +263,31 @@ def annotate_vcf( n_hom_ref_list = [] n_no_cov_list = [] any_found = False - fail_count_total = 0 for alt in variant.ALT: key = (pos, variant.REF, alt) - AC, _, in_parquet, n_fail, n_het, n_hom_alt, n_hom_ref, n_no_cov = ann.get( + AC, _, in_parquet, _, n_het, n_hom_alt, n_hom_ref, n_no_cov = ann.get( key, (0, AN, False, 0, 0, 0, 0, 0) ) if in_parquet: any_found = True ac_list.append(AC) - af_list.append(AC / AN) + af_list.append(AC / AN if AN > 0 else None) n_het_list.append(n_het) n_hom_alt_list.append(n_hom_alt) n_hom_ref_list.append(n_hom_ref) n_no_cov_list.append(n_no_cov) - fail_count_total += n_fail # cyvcf2 Number=A fields must be set as comma-separated strings variant.INFO["AFQUERY_AC"] = ",".join(str(v) for v in ac_list) variant.INFO["AFQUERY_AN"] = AN - variant.INFO["AFQUERY_AF"] = ",".join("%.10g" % v for v in af_list) + variant.INFO["AFQUERY_AF"] = ",".join( + "." if v is None else "%.10g" % v for v in af_list + ) variant.INFO["AFQUERY_N_HET"] = ",".join(str(v) for v in n_het_list) variant.INFO["AFQUERY_N_HOM_ALT"] = ",".join(str(v) for v in n_hom_alt_list) variant.INFO["AFQUERY_N_HOM_REF"] = ",".join(str(v) for v in n_hom_ref_list) - variant.INFO["AFQUERY_N_FAIL"] = fail_count_total + variant.INFO["AFQUERY_N_FAIL"] = N_FAIL variant.INFO["AFQUERY_N_NO_COVERAGE"] = ",".join(str(v) for v in n_no_cov_list) if any_found: stats["n_annotated"] += 1 diff --git a/src/afquery/dump.py b/src/afquery/dump.py index 2d062fe..4efa0ae 100644 --- a/src/afquery/dump.py +++ b/src/afquery/dump.py @@ -17,6 +17,11 @@ BUCKET_SIZE = storage.BUCKET_SIZE +def _af(stats): + """AF of one variant's stats, or None (an empty cell) when AN is 0.""" + return stats.AC / stats.AN if stats.AN > 0 else None + + def _build_groups(engine, base_sf, by_sex, by_tech, by_phenotype, all_groups): """Build list of (label, SampleFilter) for desagregation groups. @@ -186,12 +191,12 @@ def _dump_bucket_worker( for row, bitmaps in zip(rows, unpacked): pos, ref, alt = row[0], row[1], row[2] - # Base eligible / AN + # Base eligible set; a position with no eligible sample is skipped if pos not in pos_cache: pos_cache[pos] = engine._compute_eligible(chrom, pos, base_bm) - eligible, AN = pos_cache[pos] + eligible, eligible_AN = pos_cache[pos] - if AN == 0: + if eligible_AN == 0: continue s = engine._variant_stats( @@ -208,8 +213,8 @@ def _dump_bucket_worker( "ref": ref, "alt": alt, "AC": s.AC, - "AN": AN, - "AF": s.AC / AN, + "AN": s.AN, + "AF": _af(s), "N_HET": s.N_HET, "N_HOM_ALT": s.N_HOM_ALT, "N_HOM_REF": s.N_HOM_REF, @@ -223,7 +228,7 @@ def _dump_bucket_worker( cache_key = (g_idx, pos) if cache_key not in group_pos_cache: group_pos_cache[cache_key] = engine._compute_eligible(chrom, pos, g_bm) - g_eligible, g_AN = group_pos_cache[cache_key] + g_eligible = group_pos_cache[cache_key][0] g = engine._variant_stats( chrom, pos, g_eligible, bitmaps, sites[pos], @@ -231,8 +236,8 @@ def _dump_bucket_worker( ) out_row[f"AC_{label}"] = g.AC - out_row[f"AN_{label}"] = g_AN - out_row[f"AF_{label}"] = g.AC / g_AN if g_AN > 0 else 0.0 + out_row[f"AN_{label}"] = g.AN + out_row[f"AF_{label}"] = _af(g) out_row[f"N_HET_{label}"] = g.N_HET out_row[f"N_HOM_ALT_{label}"] = g.N_HOM_ALT out_row[f"N_HOM_REF_{label}"] = g.N_HOM_REF diff --git a/src/afquery/query.py b/src/afquery/query.py index 1e7c311..6c64349 100644 --- a/src/afquery/query.py +++ b/src/afquery/query.py @@ -21,11 +21,13 @@ class SiteEvidence(NamedTuple): """Carriers of any allele at one position, pooled across its Parquet rows.""" pass_bm: BitMap carrier_bm: BitMap + fail_bm: BitMap quality_pass_bm: "BitMap | None" class VariantStats(NamedTuple): AC: int + AN: int N_HET: int N_HOM_ALT: int N_HOM_REF: int @@ -323,13 +325,15 @@ def _site_evidence(unpacked_rows: list[tuple]) -> SiteEvidence: """Pool the unpacked bitmaps of every allele stored at one position.""" pass_bm = BitMap() carrier_bm = BitMap() + site_fail_bm = BitMap() quality_pass_bm = None for het_bm, hom_bm, fail_bm, _filtered_bm, qp_bm in unpacked_rows: pass_bm |= het_bm | hom_bm carrier_bm |= het_bm | hom_bm | fail_bm + site_fail_bm |= fail_bm if qp_bm is not None: quality_pass_bm = qp_bm if quality_pass_bm is None else quality_pass_bm | qp_bm - return SiteEvidence(pass_bm, carrier_bm, quality_pass_bm) + return SiteEvidence(pass_bm, carrier_bm, site_fail_bm, quality_pass_bm) @classmethod def _site_evidence_by_pos( @@ -358,26 +362,33 @@ def _variant_stats( min_observed: int, min_quality_evidence: int, ) -> VariantStats: - """Genotype tallies for one allele among the eligible samples. + """Genotype tallies, AC and AN for one allele among the eligible samples. ``site`` pools every allele at the position (None when this is the only - one). Eligible samples carrying only another allele are in none of the + one). A sample with a failed call for any allele at the position has no + reliable genotype there: it is counted in N_FAIL, shared by every allele + of the site, and left out of AC and AN. + + Eligible samples carrying only another allele are in none of the returned categories, so at a multi-allelic site N_HET + N_HOM_ALT + N_HOM_REF + N_FAIL + N_NO_COVERAGE falls short of n_eligible by exactly that number of samples. """ het_bm, hom_bm, fail_bm, filtered_bm, quality_pass_bm = bitmaps + site_fail = (site.fail_bm if site is not None else fail_bm) & eligible + called = eligible - site_fail haploid_elig, diploid_elig = split_ploidy( - eligible, self._male_bm, self._female_bm, chrom, pos, self._genome_build + called, self._male_bm, self._female_bm, chrom, pos, self._genome_build ) - het_elig = het_bm & eligible - hom_elig = hom_bm & eligible + het_elig = het_bm & called + hom_elig = hom_bm & called AC = (len((het_elig | hom_elig) & haploid_elig) + len(het_elig & diploid_elig) + 2 * len(hom_elig & diploid_elig)) + AN = len(haploid_elig) + 2 * len(diploid_elig) N_HET = len(het_elig & diploid_elig) N_HOM_ALT = len(hom_elig & diploid_elig) + len((het_elig | hom_elig) & haploid_elig) - N_FAIL = len(fail_bm & eligible) + N_FAIL = len(site_fail) no_cov_bm = self._compute_no_coverage_bm( eligible, het_bm, hom_bm, fail_bm, min_pass, min_observed, @@ -389,18 +400,18 @@ def _variant_stats( N_NO_COVERAGE = len(no_cov_bm) N_OTHER = 0 if site is not None: - N_OTHER = len((site.carrier_bm - (het_bm | hom_bm | fail_bm)) & eligible) + N_OTHER = len((site.carrier_bm - (het_bm | hom_bm | site.fail_bm)) & eligible) N_HOM_REF = len(eligible) - N_HET - N_HOM_ALT - N_FAIL - N_NO_COVERAGE - N_OTHER - return VariantStats(AC, N_HET, N_HOM_ALT, N_HOM_REF, N_FAIL, N_NO_COVERAGE) + return VariantStats(AC, AN, N_HET, N_HOM_ALT, N_HOM_REF, N_FAIL, N_NO_COVERAGE) @staticmethod def _make_result( chrom: str, pos: int, ref: str, alt: str, - AN: int, eligible: BitMap, stats: VariantStats, + eligible: BitMap, stats: VariantStats, ) -> QueryResult: return QueryResult( variant=VariantKey(chrom=chrom, pos=pos, ref=ref, alt=alt), - AC=stats.AC, AN=AN, AF=stats.AC / AN if AN > 0 else None, + AC=stats.AC, AN=stats.AN, AF=stats.AC / stats.AN if stats.AN > 0 else None, n_samples_eligible=len(eligible), N_HET=stats.N_HET, N_HOM_ALT=stats.N_HOM_ALT, N_HOM_REF=stats.N_HOM_REF, N_FAIL=stats.N_FAIL, N_NO_COVERAGE=stats.N_NO_COVERAGE, @@ -412,7 +423,12 @@ def _compute_eligible( pos: int, sample_bitmap: BitMap, ) -> tuple[BitMap, int]: - """Return (eligible_bitmap, AN) for a single position.""" + """Return (eligible_bitmap, AN of every eligible sample) for one position. + + An AN of 0 marks a position where no selected sample can be genotyped. + The AN reported for a variant comes from _variant_stats, which leaves + out samples with a failed call at the position. + """ covered = BitMap(self._covered_all_bm) # start with WGS samples (always covered) for tech_id, capture_idx in self._targeted_capture.items(): if capture_idx.covers(chrom, pos): @@ -498,7 +514,7 @@ def query(self, params: QueryParams) -> list[QueryResult]: chrom, pos, eligible, bitmaps, site, sf.min_pass, sf.min_observed, sf.min_quality_evidence, ) - results.append(self._make_result(chrom, pos, ref, alt, AN, eligible, stats)) + results.append(self._make_result(chrom, pos, ref, alt, eligible, stats)) if params.ref is not None: results = [r for r in results if r.variant.ref == params.ref] if params.alt is not None: @@ -675,12 +691,12 @@ def _query_region_inner( results = [] for row, bitmaps in zip(rows, unpacked): pos, ref, alt = row[0], row[1], row[2] - eligible, AN = pos_data[pos] + eligible = pos_data[pos][0] stats = self._variant_stats( chrom, pos, eligible, bitmaps, sites[pos], min_pass, min_observed, min_quality_evidence, ) - results.append(self._make_result(chrom, pos, ref, alt, AN, eligible, stats)) + results.append(self._make_result(chrom, pos, ref, alt, eligible, stats)) results.sort(key=lambda r: (r.variant.pos, r.variant.alt)) return results @@ -748,12 +764,12 @@ def _query_batch_inner( pos, ref, alt = row[0], row[1], row[2] if (pos, ref, alt) not in requested_variants: continue - eligible, AN = pos_data[pos] + eligible = pos_data[pos][0] stats = self._variant_stats( chrom, pos, eligible, bitmaps, sites[pos], min_pass, min_observed, min_quality_evidence, ) - results.append(self._make_result(chrom, pos, ref, alt, AN, eligible, stats)) + results.append(self._make_result(chrom, pos, ref, alt, eligible, stats)) results.sort(key=lambda r: (r.variant.pos, r.variant.alt)) return results @@ -853,9 +869,11 @@ def variant_info(self, params: QueryParams) -> list[SampleCarrier]: for row in rows: row_pos, ref, alt = row[0], row[1], row[2] het_bm, hom_bm, fail_bm, filtered_bm, quality_pass_bm = self._unpack_bitmaps(row[3:]) - het_elig = het_bm & eligible - hom_elig = hom_bm & eligible - fail_elig = fail_bm & eligible + # As in _variant_stats, a failed call for any allele at the position + # makes the sample's genotype there unknown, for every allele. + fail_elig = (site.fail_bm if site is not None else fail_bm) & eligible + het_elig = (het_bm & eligible) - fail_elig + hom_elig = (hom_bm & eligible) - fail_elig no_cov_bm = self._compute_no_coverage_bm( eligible, het_bm, hom_bm, fail_bm, sf.min_pass, sf.min_observed, @@ -869,17 +887,14 @@ def variant_info(self, params: QueryParams) -> list[SampleCarrier]: seen: set[int] = set() for sid in sorted(hom_elig): seen.add(sid) - is_fail = sid in fail_elig - carriers.append(self._make_carrier(sid, "alt" if is_fail else "hom", not is_fail)) + carriers.append(self._make_carrier(sid, "hom", True)) for sid in sorted(het_elig): if sid not in seen: seen.add(sid) - is_fail = sid in fail_elig - carriers.append(self._make_carrier(sid, "alt" if is_fail else "het", not is_fail)) + carriers.append(self._make_carrier(sid, "het", True)) for sid in sorted(fail_elig): - if sid not in seen: - seen.add(sid) - carriers.append(self._make_carrier(sid, "alt", False)) + seen.add(sid) + carriers.append(self._make_carrier(sid, "alt", False)) for sid in sorted(no_cov_elig): if sid not in seen: seen.add(sid) diff --git a/tests/conftest.py b/tests/conftest.py index dcdd57d..fab429b 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -42,14 +42,15 @@ (1, "G20"), (2, "G20"), (3, "G20"), (5, "G20"), ] -# (chrom, pos, ref, alt, het_ids, hom_ids, fail_ids) +# (chrom, pos, ref, alt, het_ids, hom_ids, fail_ids); as in a real build, a sample +# whose call failed FILTER is in fail_ids only VARIANTS = [ - ("chr1", 1500, "A", "T", [0, 5], [2], [0]), # S00 fails FILTER + ("chr1", 1500, "A", "T", [5], [2], [0]), # S00 fails FILTER ("chr1", 3500, "G", "C", [7], [], []), # no failures - ("chr1", 5000, "T", "G", [0, 1], [3], [1, 3]), # S01, S03 fail FILTER + ("chr1", 5000, "T", "G", [0], [], [1, 3]), # S01, S03 fail FILTER ("chrX", 5000000, "A", "G", [0, 2], [], []), # no failures ("chrY", 500000, "T", "C", [0, 1], [4], []), # no failures - ("chrM", 100, "C", "A", [0, 2, 5], [], [5]), # S05 fails FILTER + ("chrM", 100, "C", "A", [0, 2], [], [5]), # S05 fails FILTER ] diff --git a/tests/data/expected_results.json b/tests/data/expected_results.json index c9543cd..0c6d1f6 100644 --- a/tests/data/expected_results.json +++ b/tests/data/expected_results.json @@ -1,6 +1,6 @@ [ { - "description": "chr1:1500 E11.9 both sexes", + "description": "chr1:1500 E11.9 both sexes (S00's failed call is out of AC and AN)", "chrom": "chr1", "pos": 1500, "ref": "A", @@ -9,9 +9,9 @@ "E11.9" ], "sex": "both", - "AC": 4, - "AN": 10, - "AF": 0.4, + "AC": 3, + "AN": 8, + "AF": 0.375, "n_eligible": 5 }, { diff --git a/tests/data/vcfs/S00.vcf b/tests/data/vcfs/S00.vcf index afdd894..b04c249 100644 --- a/tests/data/vcfs/S00.vcf +++ b/tests/data/vcfs/S00.vcf @@ -1,12 +1,13 @@ ##fileformat=VCFv4.2 ##FILTER= +##FILTER= ##contig= ##contig= ##contig= ##contig= ##FORMAT= #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT S00 -chr1 1500 . A T . PASS . GT 0/1 +chr1 1500 . A T . LowQual . GT 0/1 chr1 5000 . T G . PASS . GT 0/1 chrX 5000000 . A G . PASS . GT 0/1 chrY 500000 . T C . PASS . GT 0/1 diff --git a/tests/data/vcfs/S01.vcf b/tests/data/vcfs/S01.vcf index f5fe6e5..6af6826 100644 --- a/tests/data/vcfs/S01.vcf +++ b/tests/data/vcfs/S01.vcf @@ -1,8 +1,9 @@ ##fileformat=VCFv4.2 ##FILTER= +##FILTER= ##contig= ##contig= ##FORMAT= #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT S01 -chr1 5000 . T G . PASS . GT 0/1 +chr1 5000 . T G . LowQual . GT 0/1 chrY 500000 . T C . PASS . GT 0/1 diff --git a/tests/data/vcfs/S03.vcf b/tests/data/vcfs/S03.vcf index aa94b9e..cecb57a 100644 --- a/tests/data/vcfs/S03.vcf +++ b/tests/data/vcfs/S03.vcf @@ -1,6 +1,7 @@ ##fileformat=VCFv4.2 ##FILTER= +##FILTER= ##contig= ##FORMAT= #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT S03 -chr1 5000 . T G . PASS . GT 1/1 +chr1 5000 . T G . LowQual . GT 1/1 diff --git a/tests/data/vcfs/S05.vcf b/tests/data/vcfs/S05.vcf index 3eff5ad..a83925d 100644 --- a/tests/data/vcfs/S05.vcf +++ b/tests/data/vcfs/S05.vcf @@ -1,8 +1,9 @@ ##fileformat=VCFv4.2 ##FILTER= +##FILTER= ##contig= ##contig= ##FORMAT= #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT S05 chr1 1500 . A T . PASS . GT 0/1 -chrM 100 . C A . PASS . GT 0/1 +chrM 100 . C A . LowQual . GT 0/1 diff --git a/tests/oracle.py b/tests/oracle.py index e77d740..7e8afd6 100644 --- a/tests/oracle.py +++ b/tests/oracle.py @@ -137,18 +137,30 @@ def eligible(self, chrom: str, pos: int, samples=None) -> list[str]: def expect(self, chrom, pos, ref, alt, samples=None) -> dict: """Expected query counts for one variant.""" elig = self.eligible(chrom, pos, samples) - calls = self.calls.get((chrom, pos, ref, alt), {}) + site = { + (r, a): site_calls + for (c, p, r, a), site_calls in self.calls.items() + if c == chrom and p == pos + } + calls = site.get((ref, alt), {}) + # A failed call for any allele at the position leaves the sample's + # genotype there unknown: it is out of AC and AN for every allele. + site_fail = { + s for site_calls in site.values() + for s, state in site_calls.items() if state == "fail" + } AN = 0 AC = 0 n_het = n_hom_alt = n_fail = 0 for s in elig: + if s in site_fail: + n_fail += 1 + continue haploid = self.is_haploid(s, chrom, pos) AN += 1 if haploid else 2 state = calls.get(s) - if state == "fail": - n_fail += 1 - elif state == "het": + if state == "het": if haploid: AC += 1 n_hom_alt += 1 @@ -167,10 +179,9 @@ def expect(self, chrom, pos, ref, alt, samples=None) -> dict: # n_eligible. n_other = sum( 1 for s in elig - if s not in calls and any( + if s not in site_fail and s not in calls and any( s in other_calls - for (c, p, _r, a), other_calls in self.calls.items() - if c == chrom and p == pos and (_r, a) != (ref, alt) + for key, other_calls in site.items() if key != (ref, alt) ) ) @@ -181,6 +192,7 @@ def expect(self, chrom, pos, ref, alt, samples=None) -> dict: "N_HOM_ALT": n_hom_alt, "N_FAIL": n_fail, "N_HOM_REF": len(elig) - n_het - n_hom_alt - n_fail - n_other, + "n_eligible": len(elig), } def variants(self) -> list[tuple[str, int, str, str]]: diff --git a/tests/test_annotate.py b/tests/test_annotate.py index 5241f66..872bd98 100644 --- a/tests/test_annotate.py +++ b/tests/test_annotate.py @@ -2,11 +2,14 @@ Expected values for phenotype=["E11.9"], sex="both", GRCh37: - chr1:1500 A>T — in DB — eligible={0,1,2,5,6} AN=10 AC=4 AF=0.4 + chr1:1500 A>T — in DB — eligible={0,1,2,5,6} AN=8 AC=3 AF=0.375 chr1:2000 A>G — absent — eligible={0,1,2,5,6} AN=10 AC=0 AF=0.0 chr1:3500 G>C,T — C in DB, T absent — eligible={0,1,2,7} AN=8 AC=[1,0] AF=[0.125,0.0] - chr1:5000 T>G — in DB — eligible={0,1,2} AN=6 AC=2 AF=0.333… + chr1:5000 T>G — in DB — eligible={0,1,2} AN=4 AC=1 AF=0.25 + +Samples with a failed call at a position (S00 at 1500, S01 at 5000) are left out +of AC and AN there. Coverage derivation: WGS (tech 0, samples 0-3): covers all positions @@ -59,9 +62,10 @@ def test_annotate_known_variant_1500(test_db, tmp_path): records = _read_annotated(out) r = records[1500] - assert r["AN"] == 10 - assert tuple(r["AC"]) == (4,) - assert abs(r["AF"][0] - 0.4) < 1e-6 + assert r["AN"] == 8 + assert tuple(r["AC"]) == (3,) + assert abs(r["AF"][0] - 0.375) < 1e-6 + assert r["N_FAIL"] == 1 def test_annotate_known_variant_5000(test_db, tmp_path): @@ -72,9 +76,10 @@ def test_annotate_known_variant_5000(test_db, tmp_path): records = _read_annotated(out) r = records[5000] - assert r["AN"] == 6 - assert tuple(r["AC"]) == (2,) - assert abs(r["AF"][0] - 2 / 6) < 1e-6 + assert r["AN"] == 4 + assert tuple(r["AC"]) == (1,) + assert abs(r["AF"][0] - 0.25) < 1e-6 + assert r["N_FAIL"] == 1 # --------------------------------------------------------------------------- @@ -294,8 +299,8 @@ def test_annotate_single_chrom_multi_bucket_parallel(test_db, tmp_path): parallel = _read_annotated(out_parallel) assert serial == parallel # chr1:1500 (bucket 0) is in the DB → annotated - assert serial[1500]["AN"] == 10 - assert tuple(serial[1500]["AC"]) == (4,) + assert serial[1500]["AN"] == 8 + assert tuple(serial[1500]["AC"]) == (3,) # chr1:1500000 (bucket 1) is covered by WGS but not in Parquet → AC=0, AN>0 assert serial[1500000]["AN"] == 6 assert tuple(serial[1500000]["AC"]) == (0,) diff --git a/tests/test_cli.py b/tests/test_cli.py index b89e982..8afdb9d 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -47,9 +47,9 @@ def test_query_text_format(runner, test_db): "--locus", "1:1500", "--phenotype", "E11.9", ]) assert result.exit_code == 0 - assert "AC=4" in result.output - assert "AN=10" in result.output - assert "AF=0.4000" in result.output + assert "AC=3" in result.output + assert "AN=8" in result.output + assert "AF=0.3750" in result.output def test_query_text_no_results(runner, test_db): @@ -72,9 +72,9 @@ def test_query_json_format(runner, test_db): assert result.exit_code == 0 data = json.loads(result.output) assert len(data) == 1 - assert data[0]["AC"] == 4 - assert data[0]["AN"] == 10 - assert abs(data[0]["AF"] - 0.4) < 1e-9 + assert data[0]["AC"] == 3 + assert data[0]["AN"] == 8 + assert abs(data[0]["AF"] - 0.375) < 1e-9 def test_query_json_empty(runner, test_db): @@ -100,8 +100,8 @@ def test_query_tsv_format(runner, test_db): assert lines[0].startswith("chrom\tpos\tref\talt") assert len(lines) == 2 # header + 1 variant fields = lines[1].split("\t") - assert fields[4] == "4" # AC - assert fields[5] == "10" # AN + assert fields[4] == "3" # AC + assert fields[5] == "8" # AN # --- sex filter --- diff --git a/tests/test_dump.py b/tests/test_dump.py index fc8ac78..58eae39 100644 --- a/tests/test_dump.py +++ b/tests/test_dump.py @@ -24,6 +24,14 @@ def _parse_csv(text: str) -> list[dict]: return list(reader) +def _check_ac_zero_af(row: dict, suffix: str) -> None: + """An AC=0 column set has AF=0, or an empty AF when AN=0 (only failed calls).""" + if row[f"AN{suffix}"] == "0": + assert row[f"AF{suffix}"] == "" + else: + assert float(row[f"AF{suffix}"]) == 0.0 + + def _dump_to_string(db_path, **kwargs) -> tuple[str, dict]: """Run Database.dump() into a StringIO and return (csv_text, stats).""" buf = io.StringIO() @@ -152,7 +160,8 @@ def test_ac_positive_filter(self, test_db): assert int(row["AC"]) > 0, f"AC=0 in row: {row}" def test_include_ac_zero_flag(self, test_db): - """include_ac_zero=True produces a superset; AC=0 rows have AF=0 and AN>0. + """include_ac_zero=True produces a superset; AC=0 rows have AF=0, or an + empty AF when AN=0 because every eligible sample has a failed call. chr1:3500 G>C has only S07 (phenotype J45) as carrier. When filtered by phenotype I10, no I10 sample carries this variant but several are covered @@ -174,8 +183,13 @@ def test_include_ac_zero_flag(self, test_db): ac_zero_rows = [r for r in rows_all if int(r["AC"]) == 0] assert len(ac_zero_rows) > 0, "Expected at least one AC=0 row" for row in ac_zero_rows: - assert float(row["AF"]) == 0.0 - assert int(row["AN"]) > 0 + _check_ac_zero_af(row, "") + + # chr1:5000: both eligible I10 samples (S01, S03) have a failed call there + [failed_only] = [r for r in rows_all if r["pos"] == "5000"] + assert failed_only["AN"] == "0" + assert failed_only["AF"] == "" + assert failed_only["N_FAIL"] == "2" def test_include_ac_zero_cli(self, test_db, tmp_path): """--all-variants CLI flag includes AC=0 rows for a phenotype-filtered dump.""" @@ -214,8 +228,9 @@ def test_include_ac_zero_with_disaggregation(self, test_db): ac_zero_rows = [r for r in rows if int(r["AC"]) == 0] assert len(ac_zero_rows) > 0, "Expected AC=0 rows with disaggregation" for row in ac_zero_rows: - assert float(row["AF"]) == 0.0 - assert int(row["AN"]) > 0 + _check_ac_zero_af(row, "") + for sex in ("male", "female"): + _check_ac_zero_af(row, f"_{sex}") sex_ac_cols = [k for k in row if k.startswith("AC_")] assert sex_ac_cols, "Expected disaggregated AC_ columns" for col in sex_ac_cols: diff --git a/tests/test_haploid_stats.py b/tests/test_haploid_stats.py index 0cf935a..87f3338 100644 --- a/tests/test_haploid_stats.py +++ b/tests/test_haploid_stats.py @@ -32,7 +32,7 @@ def test_chrY_all_haploid(self, test_db): assert r.AC == 2 def test_chrM_all_haploid(self, test_db): - """chrM:100 het_ids=[0,2,5] — all haploid → N_HET=0.""" + """chrM:100 het_ids=[0,2] — all haploid → N_HET=0.""" db = Database(test_db) results = db.query(chrom="chrM", pos=100) assert len(results) == 1 @@ -49,11 +49,13 @@ def test_autosome_diploid_regression(self, test_db): assert len(results) == 1 r = results[0] # Eligible: WGS{0,1,2,3} ∪ WES_A{4,5,6} (chr1:1500 covered by both) - # het_ids=[0,5], hom_ids=[2], fail_ids=[0] - # All diploid → N_HET = 2 (S00,S05), N_HOM_ALT = 1 (S02) - assert r.N_HET == 2, f"Expected N_HET=2, got {r.N_HET}" + # het_ids=[5], hom_ids=[2], fail_ids=[0] + # All diploid → N_HET = 1 (S05), N_HOM_ALT = 1 (S02), N_FAIL = 1 (S00) + assert r.N_HET == 1, f"Expected N_HET=1, got {r.N_HET}" assert r.N_HOM_ALT == 1, f"Expected N_HOM_ALT=1, got {r.N_HOM_ALT}" - assert r.AC == 4 # 2 het + 2*1 hom + assert r.N_FAIL == 1 + assert r.AC == 3 # 1 het + 2*1 hom + assert r.AN == 12 # 6 called diploid samples of the 7 eligible class TestHaploidInvariants: diff --git a/tests/test_multiallelic.py b/tests/test_multiallelic.py index 369b021..71c0463 100644 --- a/tests/test_multiallelic.py +++ b/tests/test_multiallelic.py @@ -200,3 +200,86 @@ def test_annotate_absent_allele_shares_site_coverage(panel_db, tmp_path): [v] = list(cyvcf2.VCF(str(out))) assert tuple(v.INFO.get("AFQUERY_N_NO_COVERAGE")) == (2, 2, 2) assert tuple(v.INFO.get("AFQUERY_N_HOM_REF")) == (2, 2, 2) + + +# --------------------------------------------------------------------------- +# Failed calls are judged per position +# --------------------------------------------------------------------------- + +FAIL_COHORT = [ + ("AT", "female", "WGS", [("chr1", 5000, "G", "A,T", "1/2")]), + ("AHET", "female", "WGS", [("chr1", 5000, "G", "A", "0/1")]), + ("THOM", "male", "WGS", [("chr1", 5000, "G", "T", "1/1")]), + # As written by normalize_vcf.sh: split records, GT masked at a failed site + ("F_SPLIT", "male", "WGS", [("chr1", 5000, "G", "A", "./.", "LowQual"), + ("chr1", 5000, "G", "T", "./.", "LowQual")]), + ("F_MULTI", "male", "WGS", [("chr1", 5000, "G", "A,T", "1/2", "LowQual")]), + ("F_A", "female", "WGS", [("chr1", 5000, "G", "A", "0/1", "LowQual")]), + ("REF1", "female", "WGS", [("chr1", 100, "C", "G", "0/1")]), + ("REF2", "male", "WGS", [("chr1", 100, "C", "G", "0/1")]), +] + + +@pytest.fixture(scope="module") +def fail_db(tmp_path_factory): + return _build_db(tmp_path_factory.mktemp("multiallelic_fail"), FAIL_COHORT) + + +def test_failed_calls_leave_an_for_every_allele(fail_db): + """A failed call for any allele is out of AN: AN = AC_ref + sum of AC.""" + by_alt = {r.variant.alt: r for r in fail_db.query(*SITE)} + assert (by_alt["A"].AC, by_alt["T"].AC) == (2, 3) + for r in by_alt.values(): + assert r.AN == 10 # 5 called diploid samples + assert r.N_FAIL == 3 # F_SPLIT, F_MULTI, F_A: the same for both alleles + assert r.n_samples_eligible == 8 + assert (by_alt["A"].N_HOM_REF, by_alt["T"].N_HOM_REF) == (2, 2) + ac_ref = 1 + 2 + 2 # AHET, REF1, REF2 + assert by_alt["A"].AN == ac_ref + by_alt["A"].AC + by_alt["T"].AC + + +def test_annotate_counts_each_failed_sample_once(fail_db, tmp_path): + vcf_in = tmp_path / "in.vcf" + vcf_in.write_text( + "##fileformat=VCFv4.2\n" + "##contig=\n" + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n" + "chr1\t5000\t.\tG\tA,T\t.\tPASS\t.\n" + ) + out = tmp_path / "out.vcf" + fail_db.annotate_vcf(str(vcf_in), str(out)) + [v] = list(cyvcf2.VCF(str(out))) + assert v.INFO.get("AFQUERY_N_FAIL") == 3 + assert v.INFO.get("AFQUERY_AN") == 10 + + +def test_variant_info_lists_site_failures_for_each_allele(fail_db): + """F_A failed at the position, so it is listed as alt/FAIL for T as well, + matching N_FAIL.""" + with warnings.catch_warnings(): + warnings.simplefilter("ignore", AfqueryWarning) + carriers = fail_db.variant_info(*SITE, alt="T") + failed = sorted(c.sample_name for c in carriers if c.genotype == "alt") + assert failed == ["F_A", "F_MULTI", "F_SPLIT"] + + +def test_all_eligible_failed_gives_an_zero(tmp_path): + db = _build_db(tmp_path, [ + ("F1", "female", "WGS", [("chr1", 5000, "G", "A", "0/1", "LowQual")]), + ("F2", "male", "WGS", [("chr1", 5000, "G", "A", "./.", "LowQual")]), + ]) + [r] = db.query(*SITE) + assert (r.AC, r.AN, r.AF, r.N_FAIL) == (0, 0, None, 2) + vcf_in = tmp_path / "in.vcf" + vcf_in.write_text( + "##fileformat=VCFv4.2\n" + "##contig=\n" + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\n" + "chr1\t5000\t.\tG\tA\t.\tPASS\t.\n" + ) + out = tmp_path / "out.vcf" + db.annotate_vcf(str(vcf_in), str(out)) + [v] = list(cyvcf2.VCF(str(out))) + assert v.INFO.get("AFQUERY_AN") == 0 + assert v.INFO.get("AFQUERY_N_FAIL") == 2 + assert "AFQUERY_AF=." in str(v) diff --git a/tests/test_no_coverage.py b/tests/test_no_coverage.py index a943a91..4b6565f 100644 --- a/tests/test_no_coverage.py +++ b/tests/test_no_coverage.py @@ -86,9 +86,8 @@ def test_invariant_holds_with_filter(self, test_db): class TestPhase1MinObserved: def test_min_observed_counts_fail(self, test_db): """min_observed counts fail samples as evidence — fail-only WES tech survives.""" - # chr1:1500 → WES_A has S05 het + S00 fail (S00 is WGS, not WES_A). Need careful look. # samples 0-3 are WGS, 4-6 are WES_A, 7-9 are WES_B. - # chr1:1500 het=[0,5], hom=[2], fail=[0]. WES_A carriers: het={5} (1 PASS carrier) + # chr1:1500 het=[5], hom=[2], fail=[0]. WES_A carriers: het={5} (1 PASS carrier) # WES_A's only carrier (S05) is in het → also counts as observed (1). # So min_observed=2 should filter. db = _db(test_db) diff --git a/tests/test_oracle_consistency.py b/tests/test_oracle_consistency.py index d7b97ab..a9895cc 100644 --- a/tests/test_oracle_consistency.py +++ b/tests/test_oracle_consistency.py @@ -28,7 +28,7 @@ def _assert_matches(db_dir, cohort, samples=None, phenotype=None, label=""): got = [r for r in db.query(**kwargs) if (r.variant.ref, r.variant.alt) == (ref, alt)] - if want["AN"] == 0: + if want["n_eligible"] == 0: assert not got or got[0].AC == 0 continue @@ -37,6 +37,7 @@ def _assert_matches(db_dir, cohort, samples=None, phenotype=None, label=""): actual = { "AC": r.AC, "AN": r.AN, "N_HET": r.N_HET, "N_HOM_ALT": r.N_HOM_ALT, "N_FAIL": r.N_FAIL, "N_HOM_REF": r.N_HOM_REF, + "n_eligible": r.n_samples_eligible, } assert actual == want, f"{label} {chrom}:{pos} {ref}>{alt}" checked += 1 @@ -178,6 +179,9 @@ def test_multiallelic_cohort_matches_oracle(tmp_path): "T_HET": ("female", "WGS", [("chr1", 5000, "G", "T", "0/1", "PASS")]), "AT_HET": ("male", "WGS", [("chr1", 5000, "G", "A,T", "1/2", "PASS")]), "T_FAIL": ("female", "WGS", [("chr1", 5000, "G", "T", "0/1", "LowQual")]), + # A failed multi-allelic site as normalize_vcf.sh writes it + "F_SPLIT": ("male", "WGS", [("chr1", 5000, "G", "A", "./.", "LowQual"), + ("chr1", 5000, "G", "T", "./.", "LowQual")]), "REF": ("female", "WGS", [("chr1", 100, "C", "G", "0/1", "PASS")]), "P_T": ("female", "PANEL", [("chr1", 5000, "G", "T", "1/1", "PASS")]), "P_REF": ("male", "PANEL", [("chr1", 100, "C", "G", "0/1", "PASS")]), @@ -213,11 +217,13 @@ def test_multiallelic_cohort_matches_oracle(tmp_path): cohort = oracle.Cohort(manifest_path, bed_dir=beds) _assert_matches(db, cohort, label="multiallelic") - # Spot-check chr1:5000, where all 11 samples are eligible (the panel BED + # Spot-check chr1:5000, where all 12 samples are eligible (the panel BED # covers it). Hom-ref for both alleles: REF, P_REF, P_OFF, X_A, X_C. got = {r.variant.alt: r for r in Database(str(db)).query("chr1", 5000)} - assert got["A"].n_samples_eligible == 11 - # A: A_HET, AT_HET het; A_HOM hom; T_HET, T_FAIL, P_T carry only T. - assert (got["A"].N_HET, got["A"].N_HOM_ALT, got["A"].N_HOM_REF) == (2, 1, 5) - # T: T_HET, AT_HET het; P_T hom; T_FAIL failed; A_HET, A_HOM carry only A. - assert (got["T"].N_HET, got["T"].N_HOM_ALT, got["T"].N_FAIL, got["T"].N_HOM_REF) == (2, 1, 1, 5) + assert got["A"].n_samples_eligible == 12 + # T_FAIL and F_SPLIT failed at the position: N_FAIL for both alleles, out of AN. + # A: A_HET, AT_HET het; A_HOM hom; T_HET, P_T carry only T. + assert (got["A"].N_HET, got["A"].N_HOM_ALT, got["A"].N_FAIL, got["A"].N_HOM_REF) == (2, 1, 2, 5) + # T: T_HET, AT_HET het; P_T hom; A_HET, A_HOM carry only A. + assert (got["T"].N_HET, got["T"].N_HOM_ALT, got["T"].N_FAIL, got["T"].N_HOM_REF) == (2, 1, 2, 5) + assert got["A"].AN == got["T"].AN == 20 diff --git a/tests/test_sample_filter.py b/tests/test_sample_filter.py index bfe19ad..ad1343b 100644 --- a/tests/test_sample_filter.py +++ b/tests/test_sample_filter.py @@ -11,7 +11,7 @@ # Test DB setup (from conftest.py): # SAMPLES: {0,1,2,3}=WGS, {4,5,6}=WES_kit_A, {7,8,9}=WES_kit_B # E11.9={0,1,2,5,6,7}, I10={1,3,5,8,9}, C50={2,4,6}, J45={0,3,7,8,9}, G20={1,2,3,5} -# chr1 pos=1500: het=[0,5], hom=[2] +# chr1 pos=1500: het=[5], hom=[2], fail=[0] # BED coverage: WES_kit_A covers chr1 1000-1999, WES_kit_B covers chr1 3000-3999 # So covered at 1500: WGS={0,1,2,3} + WES_kit_A={4,5,6} = {0,1,2,3,4,5,6} @@ -62,7 +62,7 @@ def test_sample_filter_parse_empty(): def test_filter_phenotype_exclusion(test_db): # Exclude E11.9={0,1,2,5,6,7} → remaining {3,4,8,9} # covered at 1500: {0..6} → eligible = {3,4} - # het=[0,5], hom=[2] → none in {3,4} → AC=0, AN=4 (2 diploid samples) + # het=[5], hom=[2], fail=[0] → none in {3,4} → AC=0, AN=4 (2 diploid samples) db = Database(test_db) results = db.query(chrom="chr1", pos=1500, phenotype=["^E11.9"], sex="both") assert len(results) > 0, "Should return results when AN > 0" @@ -153,7 +153,8 @@ def test_filter_combined(test_db): # tech WGS → {0,1,2,3} # sample_bm = {0,2,6,7} & {0,1,2,3} = {0,2} # covered at 1500: {0..6} → eligible = {0,2} - # het=[0,5]→{0} in eligible, hom=[2]→{2} in eligible → AC = 1 + 2 = 3 + # fail=[0]→{0}: S00's call failed → out of AC and AN + # hom=[2]→{2} in eligible → AC = 2, AN = 2 db = Database(test_db) results = db.query( chrom="chr1", pos=1500, @@ -162,8 +163,9 @@ def test_filter_combined(test_db): assert len(results) > 0 r = results[0] assert r.n_samples_eligible == 2 # {0,2} - assert r.AC == 3 # sample 0 het(+1) + sample 2 hom(+2) - assert r.AN == 4 # 2 diploid samples + assert r.AC == 2 # sample 2 hom(+2) + assert r.AN == 2 # sample 2 only; sample 0 has a failed call + assert r.N_FAIL == 1 # --------------------------------------------------------------------------- diff --git a/tests/test_update.py b/tests/test_update.py index 4d3b5eb..01fd490 100644 --- a/tests/test_update.py +++ b/tests/test_update.py @@ -172,7 +172,7 @@ def test_add_samples_new_variant(fresh_db, tmp_path): def test_add_samples_existing_variant(fresh_db, tmp_path): - # chr1:1500 A>T has het=[0,5], hom=[2] in the base DB + # chr1:1500 A>T has het=[5], hom=[2], fail=[0] in the base DB # S10 (male, wgs, E11.9) also carries chr1:1500 het → AC should increase db_before = Database(fresh_db) results_before = db_before.query(chrom="chr1", pos=1500, phenotype=["E11.9"], sex="both") @@ -257,18 +257,17 @@ def test_remove_samples_id_not_reused(fresh_db, tmp_path): def test_remove_samples_query_impact(fresh_db): - # S00 (id=0, male, wgs, E11.9+J45) is het-carrier of chr1:1500 A>T + # S00 (id=0, male, wgs, E11.9+J45) has a failed call at chr1:1500 A>T db_before = Database(fresh_db) - results_before = db_before.query(chrom="chr1", pos=1500, phenotype=["E11.9"], sex="both") - ac_before = results_before[0].AC if results_before else 0 + [before] = db_before.query(chrom="chr1", pos=1500, phenotype=["E11.9"], sex="both") remove_samples(fresh_db, ["S00"]) db_after = Database(fresh_db) - results_after = db_after.query(chrom="chr1", pos=1500, phenotype=["E11.9"], sex="both") - ac_after = results_after[0].AC if results_after else 0 + [after] = db_after.query(chrom="chr1", pos=1500, phenotype=["E11.9"], sex="both") - assert ac_after < ac_before + assert after.N_FAIL == before.N_FAIL - 1 + assert after.n_samples_eligible == before.n_samples_eligible - 1 def test_remove_nonexistent_raises(fresh_db): diff --git a/tests/test_variant_info.py b/tests/test_variant_info.py index 3c71661..a09bb43 100644 --- a/tests/test_variant_info.py +++ b/tests/test_variant_info.py @@ -13,12 +13,12 @@ # Reference data from conftest.py # --------------------------------------------------------------------------- # VARIANTS = [ -# ("chr1", 1500, "A", "T", het=[0,5], hom=[2], fail=[0]), # S00 in both het+fail +# ("chr1", 1500, "A", "T", het=[5], hom=[2], fail=[0]), # ("chr1", 3500, "G", "C", het=[7], hom=[], fail=[]), -# ("chr1", 5000, "T", "G", het=[0,1], hom=[3], fail=[1,3]), +# ("chr1", 5000, "T", "G", het=[0], hom=[], fail=[1,3]), # ("chrX", 5000000, "A", "G", het=[0,2], hom=[], fail=[]), # ("chrY", 500000, "T", "C", het=[0,1], hom=[4], fail=[]), -# ("chrM", 100, "C", "A", het=[0,2,5], hom=[], fail=[5]), +# ("chrM", 100, "C", "A", het=[0,2], hom=[], fail=[5]), # ] # SAMPLES (id, name, sex, tech_id): # 0=S00/male/WGS, 1=S01/male/WGS, 2=S02/female/WGS, 3=S03/female/WGS, @@ -51,7 +51,7 @@ def test_hom_genotype(test_db): def test_fail_genotype(test_db): - """S00 is in both het_bitmap and fail_bitmap → should be reported as alt/FAIL.""" + """S00's call failed FILTER → reported as alt/FAIL.""" db = Database(test_db) carriers = db.variant_info("chr1", 1500, ref="A", alt="T") by_id = {c.sample_id: c for c in carriers} @@ -110,23 +110,22 @@ def test_phenotypes_sorted(test_db): # --------------------------------------------------------------------------- # 3. All fail_bitmap cases — chr1:5000 T>G -# het=[0,1], hom=[3], fail=[1,3] -# → S01: het+fail → alt/FAIL -# → S03: hom+fail → alt/FAIL -# → S00: het only → het/PASS +# het=[0], hom=[], fail=[1,3] +# → S01, S03: failed calls → alt/FAIL +# → S00: het → het/PASS # --------------------------------------------------------------------------- def test_fail_chr1_5000(test_db): db = Database(test_db) carriers = db.variant_info("chr1", 5000, ref="T", alt="G") by_id = {c.sample_id: c for c in carriers} - # S00: het, no fail + # S00: het assert by_id[0].genotype == "het" assert by_id[0].filter_pass is True - # S01: het + fail + # S01: failed call (0/1 in its VCF) assert by_id[1].genotype == "alt" assert by_id[1].filter_pass is False - # S03: hom + fail + # S03: failed call (1/1 in its VCF) assert by_id[3].genotype == "alt" assert by_id[3].filter_pass is False From b712cb496797fee4808f9bd41f797420a4c94c15 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20L=C3=B3pez=20L=C3=B3pez?= Date: Wed, 14 Oct 2026 12:52:09 +0200 Subject: [PATCH 04/10] feat(query): add --fail-as-ref to count failed calls as reference Failed calls are now left out of AN. Some analyses still want the previous estimate, in which every failed sample is a hom-ref non-carrier and AF is a conservative lower bound. --fail-as-ref on query, annotate and dump (fail_as_ref in the Python API and SampleFilter) counts those samples as hom-ref in AN. AC, N_FAIL and the genotype tallies are unchanged. The AFQUERY_AN header of annotated VCFs states which convention was used. --- docs/advanced/filter-pass-tracking.md | 18 +++++++++- docs/getting-started/understanding-output.md | 2 +- docs/reference/cli.md | 3 ++ docs/reference/glossary.md | 2 +- docs/reference/python-api.md | 16 +++++++++ src/afquery/annotate.py | 12 ++++--- src/afquery/cli.py | 20 ++++++----- src/afquery/database.py | 37 ++++++++++++++++---- src/afquery/dump.py | 3 ++ src/afquery/models.py | 3 ++ src/afquery/query.py | 21 +++++++---- tests/oracle.py | 11 ++++-- tests/test_cli.py | 12 +++++++ tests/test_cli_docs_consistency.py | 5 +++ tests/test_multiallelic.py | 35 ++++++++++++++++++ tests/test_oracle_consistency.py | 12 +++++-- 16 files changed, 178 insertions(+), 34 deletions(-) diff --git a/docs/advanced/filter-pass-tracking.md b/docs/advanced/filter-pass-tracking.md index 6fff8e6..93007f6 100644 --- a/docs/advanced/filter-pass-tracking.md +++ b/docs/advanced/filter-pass-tracking.md @@ -31,7 +31,7 @@ afquery query --db ./db/ --locus chr1:925952 ``` ``` -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 +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 ``` 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. @@ -42,6 +42,22 @@ If all eligible samples have a failed call at a position, AN is 0 and AF is unde --- +## Counting failed calls as reference + +`--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 query --db ./db/ --locus chr1:925952 --fail-as-ref +``` + +``` +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. + +--- + ## Finding the failed samples `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). diff --git a/docs/getting-started/understanding-output.md b/docs/getting-started/understanding-output.md index bc74475..4d31b50 100644 --- a/docs/getting-started/understanding-output.md +++ b/docs/getting-started/understanding-output.md @@ -9,7 +9,7 @@ This page explains what each field in AFQuery output means and how to interpret | Field | Type | Description | |-------|------|-------------| | **AC** | int | Allele count — number of alt allele copies in eligible samples | -| **AN** | int | Allele number — alleles of the eligible samples with a usable call at this position, adjusted for ploidy. Samples in `N_FAIL` are left out | +| **AN** | int | Allele number — alleles of the eligible samples with a usable call at this position, adjusted for ploidy. Samples in `N_FAIL` are left out (counted as hom-ref with `--fail-as-ref`) | | **AF** | float | Allele frequency — `AC / AN`. `NA` (text/TSV) or `null` (JSON) 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). | diff --git a/docs/reference/cli.md b/docs/reference/cli.md index 4f7d707..04da002 100644 --- a/docs/reference/cli.md +++ b/docs/reference/cli.md @@ -70,6 +70,7 @@ Exactly one of `--locus`, `--region`, or `--from-file` must be provided. | `--min-pass` | INTEGER | `0` | Min PASS carriers (`het\|hom`) per partially-covered tech for hom-ref to be assumed. Non-carriers move to `N_NO_COVERAGE` if a tech falls below the threshold. (0 = disabled) | | `--min-observed` | INTEGER | `0` | Min any-VCF entries (`het\|hom\|fail`) per partially-covered tech for hom-ref to be assumed. (0 = disabled) | | `--min-quality-evidence` | INTEGER | `0` | Min quality-passing carriers per partially-covered tech. Requires a database built with `--min-dp`, `--min-gq`, `--min-qual`, or `--min-covered`. (0 = disabled) | +| `--fail-as-ref` | flag | False | Count samples with a failed call (FILTER≠PASS) as hom-ref in AN instead of leaving them out. See [FILTER=PASS Tracking](../advanced/filter-pass-tracking.md#counting-failed-calls-as-reference). | | `--no-warn` | flag | False | Suppress warnings for unknown phenotypes, technologies, and chromosomes | --- @@ -95,6 +96,7 @@ afquery annotate [OPTIONS] | `--min-pass` | INTEGER | `0` | Min PASS carriers per partially-covered tech for hom-ref to be assumed (0 = disabled) | | `--min-observed` | INTEGER | `0` | Min any-VCF entries per partially-covered tech (0 = disabled) | | `--min-quality-evidence` | INTEGER | `0` | Min quality-passing carriers per partially-covered tech. Requires a database built with `--min-dp`, `--min-gq`, `--min-qual`, or `--min-covered`. (0 = disabled) | +| `--fail-as-ref` | flag | False | Count samples with a failed call (FILTER≠PASS) as hom-ref in AN instead of leaving them out. See [FILTER=PASS Tracking](../advanced/filter-pass-tracking.md#counting-failed-calls-as-reference). | | `--no-warn` | flag | False | Suppress warnings for unknown phenotypes, technologies, and chromosomes | The annotated VCF gains `AFQUERY_AC`, `AFQUERY_AN`, `AFQUERY_AF`, @@ -131,6 +133,7 @@ afquery dump [OPTIONS] | `--min-pass` | INTEGER | `0` | Min PASS carriers per partially-covered tech for hom-ref to be assumed (0 = disabled) | | `--min-observed` | INTEGER | `0` | Min any-VCF entries per partially-covered tech (0 = disabled) | | `--min-quality-evidence` | INTEGER | `0` | Min quality-passing carriers per partially-covered tech. Requires a database built with `--min-dp`, `--min-gq`, `--min-qual`, or `--min-covered`. (0 = disabled) | +| `--fail-as-ref` | flag | False | Count samples with a failed call (FILTER≠PASS) as hom-ref in AN instead of leaving them out. See [FILTER=PASS Tracking](../advanced/filter-pass-tracking.md#counting-failed-calls-as-reference). | CSV output adds an `N_NO_COVERAGE` column (and per-group variants `N_NO_COVERAGE_