From ae50d27a310e701f89f726765835d53a348463bc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20L=C3=B3pez=20L=C3=B3pez?= Date: Wed, 30 Sep 2026 12:40:38 +0200 Subject: [PATCH 1/5] fix(variant-info): restrict carriers to samples covered at the position variant-info intersected the per-variant bitmaps with the sample filter only, so a sample with a call outside its own capture BED was listed as a carrier while query() excluded it from AC, AN and the genotype tallies. Use the same capture- and ploidy-aware eligible set as query(), so both commands agree. On the CSVS hs37d5 database this accounts for every reported mismatch, e.g. at 1:6508625 variant-info listed 120 het / 8 hom / 7 alt against N_HET=95 / N_HOM_ALT=7 / N_FAIL=0 from query. --- docs/advanced/debugging-results.md | 2 +- docs/guides/variant-info.md | 19 +++--- src/afquery/query.py | 11 +-- tests/test_variant_info.py | 106 +++++++++++++++++++++++++++++ 4 files changed, 121 insertions(+), 17 deletions(-) diff --git a/docs/advanced/debugging-results.md b/docs/advanced/debugging-results.md index 22074ee..3de8794 100644 --- a/docs/advanced/debugging-results.md +++ b/docs/advanced/debugging-results.md @@ -16,7 +16,7 @@ AN=0 means no eligible samples at the queried position. Work through these check | Position exists in database | `afquery query --db ./db/ --locus chr1:12345678` | If no result at all, the variant was not observed in any sample during ingestion. | | BED coverage (WES) | `afquery info --db ./db/` | If all eligible samples are WES and the position is outside capture regions, AN=0 is correct. | | Sample filter too restrictive | Remove `--phenotype` and `--sex` filters | Query with no filters first. If AN>0 without filters, the filter is excluding all samples. | -| WES samples missing from AN | Compare `variant-info` carriers against `query` counts | If `variant-info` shows many WES carriers but AN reflects only the WGS samples, that technology's capture index is not matching the position. `afquery` warns on open when a capture BED is empty or matches no known chromosome, naming the offending contigs. | +| WES samples missing from AN | `afquery query --db ./db/ --locus chr1:12345678 --tech WES_kit_A` | If `n_eligible=0` for a technology whose kit should cover the position, that technology's capture index is not matching it. `afquery` warns on open when a capture BED is empty or matches no known chromosome, naming the offending contigs. | | AN=0 only on chrM | `afquery query --db ./db/ --locus chrM:3243` | Mitochondrial coverage is easy to lose on its own: a BED naming it `chrMT` still matches on the autosomes, so nothing else looks wrong. Names are normalized (`MT`, `M`, `chrMT` → `chrM`) since v0.4.0 — databases built before that may need the capture BED re-checked. | | Technology filter | Remove `--tech` filter | Check if any samples match the requested technology. | diff --git a/docs/guides/variant-info.md b/docs/guides/variant-info.md index ab40333..f0953ae 100644 --- a/docs/guides/variant-info.md +++ b/docs/guides/variant-info.md @@ -15,17 +15,14 @@ afquery variant-info --db ./db/ --locus chr1:925952 !!! tip `variant-info` is the natural next step after `query` — once you find a variant of interest, use it to see which specific samples carry it. -!!! note "Carrier counts may exceed the counts reported by `query`" - `variant-info` lists every carrier present in the source VCFs. `query` additionally - restricts to samples whose capture regions cover the position, so a WES sample with a - call just outside its own BED appears here but is not counted in `AN`, `AC` or the - genotype tallies. The two commands answer different questions — "who carries it?" - versus "what is the frequency among samples that could have been called?" — so a - modest difference is expected. - - A *large* gap is worth investigating: if `variant-info` shows many carriers from a - technology that contributes nothing to `AN`, check that technology's capture index - (see [Debugging Results](../advanced/debugging-results.md)). +!!! note "Carriers match the counts reported by `query`" + `variant-info` lists only samples whose capture regions cover the position, the same + eligible set `query` uses. A WES sample with a call outside its own BED is not counted + in `AN`, `AC` or the genotype tallies, and is not listed here either. On diploid + regions the number of `het`, `hom` and `alt` carriers therefore equals `N_HET`, + `N_HOM_ALT` and `N_FAIL` for the same allele and filters. (In haploid regions a + single-allele call is listed as `het` but counted in `N_HOM_ALT`.) Earlier versions + also listed off-target carriers, so the two commands could disagree. By default all samples are queried and results are printed as an aligned text table: diff --git a/src/afquery/query.py b/src/afquery/query.py index 4939c30..b2aac57 100644 --- a/src/afquery/query.py +++ b/src/afquery/query.py @@ -780,7 +780,8 @@ def variant_info(self, params: QueryParams) -> list[SampleCarrier]: AfqueryWarning, stacklevel=3, ) - # Compute eligible (BED-aware) for no_coverage assessment + # Same eligible set as query(): a call outside the sample's capture + # region is not counted there, so it is not listed here either. eligible, _AN = self._compute_eligible(chrom, pos, sample_bm) sf = params.filter @@ -788,9 +789,9 @@ 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 & sample_bm - hom_elig = hom_bm & sample_bm - fail_elig = fail_bm & sample_bm + het_elig = het_bm & eligible + hom_elig = hom_bm & eligible + fail_elig = fail_bm & eligible no_cov_bm = self._compute_no_coverage_bm( eligible, het_bm, hom_bm, fail_bm, sf.min_pass, sf.min_observed, @@ -798,7 +799,7 @@ def variant_info(self, params: QueryParams) -> list[SampleCarrier]: quality_pass_bm=quality_pass_bm, min_quality_evidence=sf.min_quality_evidence, ) - no_cov_elig = no_cov_bm & sample_bm + no_cov_elig = no_cov_bm & eligible seen: set[int] = set() for sid in sorted(hom_elig): diff --git a/tests/test_variant_info.py b/tests/test_variant_info.py index 4ff0b24..fc6bea9 100644 --- a/tests/test_variant_info.py +++ b/tests/test_variant_info.py @@ -433,3 +433,109 @@ def test_empty_result_no_warning(test_db): # Should not warn about multiple alleles when explicitly filtered warns = [x for x in w if "alleles" in str(x.message).lower()] assert len(warns) == 0 + + +# --------------------------------------------------------------------------- +# 15. Capture regions: carriers outside the sample's BED are not listed +# --------------------------------------------------------------------------- + +def test_out_of_capture_carrier_excluded_chrY(test_db): + """chrY:500000 hom=[4]: S04 is WES_kit_A, whose BED does not cover chrY.""" + db = Database(test_db) + carriers = db.variant_info("chrY", 500000, ref="T", alt="C") + assert {c.sample_id for c in carriers} == {0, 1} + + +def test_out_of_capture_carrier_excluded_chrM(test_db): + """chrM:100 het=[0,2,5]: S05 is WES_kit_A, whose BED does not cover chrM.""" + db = Database(test_db) + carriers = db.variant_info("chrM", 100, ref="C", alt="A") + assert {c.sample_id for c in carriers} == {0, 2} + + +def _write_vcf(path, sample_name, records): + contigs = sorted({r[0] for r in records}) or ["chr1"] + with open(path, "w") as f: + f.write("##fileformat=VCFv4.2\n") + f.write('##FILTER=\n') + f.write('##FILTER=\n') + for contig in contigs: + f.write(f"##contig=\n") + f.write('##FORMAT=\n') + f.write(f"#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t{sample_name}\n") + for rec in records: + chrom, pos, ref, alt, gt = rec[:5] + flt = rec[5] if len(rec) > 5 else "PASS" + 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}.""" + from afquery.preprocess import run_preprocess + + bed_dir = tmp_path / "beds" + bed_dir.mkdir() + for tech, text in (beds or {}).items(): + (bed_dir / f"{tech}.bed").write_text(text) + rows = ["sample_name\tsex\ttech_name\tvcf_path\tphenotype_codes"] + for name, sex, tech, records in sample_defs: + vcf_path = tmp_path / f"{name}.vcf" + _write_vcf(vcf_path, name, records) + rows.append(f"{name}\t{sex}\t{tech}\t{vcf_path}\tE11.9") + manifest = tmp_path / "manifest.tsv" + manifest.write_text("\n".join(rows) + "\n") + db_path = tmp_path / "db" + db_path.mkdir() + run_preprocess( + manifest_path=str(manifest), output_dir=str(db_path), + genome_build="GRCh37", threads=1, bed_dir=str(bed_dir), + ) + return Database(str(db_path)) + + +def test_off_target_call_not_listed(tmp_path): + """A panel sample with a call outside its BED is neither counted nor listed.""" + db = _build_db(tmp_path, [ + ("WGS_HET", "female", "WGS", [("chr1", 5000, "G", "A", "0/1")]), + ("PANEL_OFF", "female", "PANEL", [("chr1", 5000, "G", "A", "0/1")]), + ("PANEL_IN", "female", "PANEL", [("chr1", 500, "C", "T", "0/1")]), + ], beds={"PANEL": "chr1\t0\t1000\n"}) + [r] = db.query("chr1", 5000) + assert r.n_samples_eligible == 1 + carriers = db.variant_info("chr1", 5000, ref="G", alt="A") + assert [c.sample_name for c in carriers] == ["WGS_HET"] + + +def test_carrier_counts_match_query(tmp_path): + """variant_info het/hom/alt counts equal query N_HET/N_HOM_ALT/N_FAIL.""" + db = _build_db(tmp_path, [ + ("W1", "female", "WGS", [("chr1", 5000, "G", "A", "0/1")]), + ("W2", "male", "WGS", [("chr1", 5000, "G", "A", "1/1")]), + ("W3", "female", "WGS", [("chr1", 5000, "G", "A", "0/1", "LowQual")]), + ("P1", "female", "PANEL", [("chr1", 5000, "G", "A", "1/1")]), + ("P2", "male", "PANEL", [("chr1", 5000, "G", "A", "0/1", "LowQual")]), + ("P3", "female", "PANEL", [("chr1", 500, "G", "A", "0/1")]), + ("P4", "female", "PANEL", [("chr1", 500, "G", "A", "1/1")]), + ], beds={"PANEL": "chr1\t0\t1000\n"}) + for pos in (500, 5000): + [r] = db.query("chr1", pos) + genotypes = [c.genotype for c in db.variant_info("chr1", pos, ref="G", alt="A")] + assert genotypes.count("het") == r.N_HET + assert genotypes.count("hom") == r.N_HOM_ALT + assert genotypes.count("alt") == r.N_FAIL + + +def test_multiallelic_1_2_sample_listed_for_each_allele(tmp_path): + """A 1/2 sample carries both alleles and is listed once per allele.""" + db = _build_db(tmp_path, [ + ("S12", "female", "WGS", [("chr1", 5000, "G", "A,T", "1/2")]), + ("SA", "female", "WGS", [("chr1", 5000, "G", "A", "0/1")]), + ]) + names_a = [c.sample_name for c in db.variant_info("chr1", 5000, alt="A")] + names_t = [c.sample_name for c in db.variant_info("chr1", 5000, alt="T")] + assert names_a == ["S12", "SA"] + assert names_t == ["S12"] + with warnings.catch_warnings(): + warnings.simplefilter("ignore", AfqueryWarning) + both = [c.sample_name for c in db.variant_info("chr1", 5000)] + assert sorted(both) == ["S12", "S12", "SA"] From e75308ccb93cd10dcd8dc77301a106260e358cc6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20L=C3=B3pez=20L=C3=B3pez?= Date: Wed, 30 Sep 2026 12:54:24 +0200 Subject: [PATCH 2/5] fix(query): exclude other-allele carriers from N_HOM_REF at multi-allelic sites N_HOM_REF was the residual n_eligible - N_HET - N_HOM_ALT - N_FAIL - N_NO_COVERAGE. The eligible set depends only on the position, so at a site with several ALT alleles every sample carrying another allele was counted as hom-ref for this one. At chr13:44995309 G>A,T this reported N_HOM_REF=165 for G>A while 36 of those samples carry G>T. The tallies for one allele now leave out eligible samples that carry only another allele at the position. At a multi-allelic site the five categories therefore add up to less than n_eligible, by exactly that number of samples; biallelic sites are unchanged. AC, AN and AF are not affected. The documented invariant is updated accordingly. Coverage evidence (--min-pass, --min-observed, --min-quality-evidence) is now judged per position: a call for any allele shows the position was sequenced, and such a sample is never reported as N_NO_COVERAGE. The per-allele computation, previously repeated in query, region, batch, dump and annotate, now lives in QueryEngine._variant_stats. annotate also excludes carriers of stored alleles when the requested allele is absent from the database. The test oracle had the same residual and is corrected; a multi-allelic cohort is added to the oracle tests. --- docs/advanced/coverage-evidence.md | 8 +- docs/advanced/ploidy-and-sex-chroms.md | 4 +- docs/getting-started/understanding-output.md | 17 +- docs/guides/annotate-vcf.md | 2 +- docs/guides/query.md | 7 +- docs/reference/python-api.md | 5 +- src/afquery/annotate.py | 50 ++-- src/afquery/dump.py | 94 ++------ src/afquery/query.py | 231 +++++++++++-------- tests/oracle.py | 14 +- tests/test_multiallelic.py | 153 ++++++++++++ tests/test_oracle_consistency.py | 59 +++++ 12 files changed, 438 insertions(+), 206 deletions(-) create mode 100644 tests/test_multiallelic.py diff --git a/docs/advanced/coverage-evidence.md b/docs/advanced/coverage-evidence.md index 364c1f2..762bc54 100644 --- a/docs/advanced/coverage-evidence.md +++ b/docs/advanced/coverage-evidence.md @@ -17,12 +17,16 @@ into `N_HOM_REF`. The flags below decide *which* samples land there. ## What `N_NO_COVERAGE` represents `N_NO_COVERAGE` counts eligible samples whose hom-ref status is not trusted -under the active criteria. The genotype invariant becomes: +under the active criteria. At a biallelic site the genotype invariant becomes: ``` N_HET + N_HOM_ALT + N_HOM_REF + N_FAIL + N_NO_COVERAGE = n_eligible ``` +At a multi-allelic site, eligible samples that carry only another ALT allele at +the position fall in none of these categories for this allele, so the sum is +lower than `n_eligible` by exactly that number of samples. + Samples in `N_NO_COVERAGE` remain in `eligible` and contribute to `AN` (just like `N_FAIL`), so AC/AN/AF stay conservative — the field never inflates allele frequencies. Two rules always hold: @@ -45,7 +49,7 @@ position. They run at query time, so no database rebuild is needed. | Flag | Effect | |------|--------| -| `--min-pass K` | A partially-covered tech must have ≥K PASS carriers (`het ∪ hom`) at the position. If it falls short, all of its non-carrier samples move from `N_HOM_REF` to `N_NO_COVERAGE`. | +| `--min-pass K` | A partially-covered tech must have ≥K PASS carriers (`het ∪ hom`) at the position. Carriers of any ALT allele there count, since each shows the position was sequenced. If it falls short, all of its non-carrier samples move from `N_HOM_REF` to `N_NO_COVERAGE`. | | `--min-observed K` | Same shape, but counts every recorded carrier (`het ∪ hom ∪ fail`). Useful when a non-PASS call still proves the position was sequenced. | When both flags are >0, both must hold (AND). The default `0` disables the diff --git a/docs/advanced/ploidy-and-sex-chroms.md b/docs/advanced/ploidy-and-sex-chroms.md index 71108e1..36acbcf 100644 --- a/docs/advanced/ploidy-and-sex-chroms.md +++ b/docs/advanced/ploidy-and-sex-chroms.md @@ -103,11 +103,11 @@ afquery query --db ./db/ --locus chrM:3243 ### Counting Identity -For every query result, the following identity holds: +For every query result at a biallelic site, the following identity holds: **N_HET + N_HOM_ALT + N_HOM_REF + N_FAIL + N_NO_COVERAGE = n_eligible** -This can be used to validate results. N_HOM_REF is the number of eligible samples that are homozygous reference (i.e., do not carry the alt allele and passed quality filters). N_NO_COVERAGE is 0 unless a coverage-evidence filter is active — see [Coverage Evidence](coverage-evidence.md). +This can be used to validate results. At a multi-allelic site, eligible samples that carry only another ALT allele at the position fall in none of these categories, so the sum is lower by exactly that number of samples (see [Multi-allelic sites](../getting-started/understanding-output.md#multi-allelic-sites)). N_HOM_REF is the number of eligible samples that are homozygous reference (i.e., carry no ALT allele at the position and passed quality filters). N_NO_COVERAGE is 0 unless a coverage-evidence filter is active — see [Coverage Evidence](coverage-evidence.md). !!! note "Mutual exclusivity" N_HET, N_HOM_ALT, N_HOM_REF, N_FAIL, and N_NO_COVERAGE are mutually exclusive. A sample with a non-ref allele but FILTER≠PASS is counted in N_FAIL only — it does not appear in N_HET or N_HOM_ALT. Likewise, N_HOM_REF counts only PASS-filtered samples. diff --git a/docs/getting-started/understanding-output.md b/docs/getting-started/understanding-output.md index e1528b9..b7d7781 100644 --- a/docs/getting-started/understanding-output.md +++ b/docs/getting-started/understanding-output.md @@ -13,7 +13,7 @@ This page explains what each field in AFQuery output means and how to interpret | **AF** | float | Allele frequency — `AC / AN`. `None` when AN=0 | | **N_HET** | int | Number of eligible samples heterozygous for the alt allele (GT=0/1) | | **N_HOM_ALT** | int | Number of eligible samples homozygous for the alt allele (GT=1/1 or GT=1). Includes haploid carriers on sex chromosomes and chrM. See [Ploidy](../advanced/ploidy-and-sex-chroms.md#genotype-counting). | -| **N_HOM_REF** | int | Number of eligible samples homozygous reference (GT=0/0 or GT=0) | +| **N_HOM_REF** | int | Number of eligible samples homozygous reference (GT=0/0 or GT=0). Samples carrying another ALT allele at the same position are not counted (see [Multi-allelic sites](#multi-allelic-sites)). | | **n_eligible** | int | Number of eligible samples — those passing the sex/phenotype/tech filters *and* covered at this position | | **N_FAIL** | int | Number of eligible samples whose call at this position had FILTER≠PASS. These samples are counted *only* in N_FAIL — not in N_HET, N_HOM_ALT, or N_HOM_REF — but they stay eligible and still count toward AN. | | **N_NO_COVERAGE** | int | Number of eligible samples whose tech lacks coverage evidence at this position. Excluded from `N_HOM_REF` to keep AC/AN conservative. Always `0` unless a coverage-evidence filter is active. See [Coverage Evidence](../advanced/coverage-evidence.md). | @@ -84,6 +84,19 @@ AN=0 means no eligible samples have coverage at this position. This happens when !!! warning "AN=0 does not mean the variant is absent" AN=0 means AFQuery has no data to compute frequency. It is not evidence of rarity. +### Multi-allelic sites + +Each ALT allele at a position is reported on its own line, and its counts refer to +that allele only. A sample carrying a different ALT allele at the same position +(for example `G/T` on the `G>A` line) is neither a carrier nor homozygous +reference for this allele, so it appears in none of `N_HET`, `N_HOM_ALT`, +`N_HOM_REF`, `N_FAIL` or `N_NO_COVERAGE`. The five counts add up to +`n_eligible` at biallelic sites; at a multi-allelic site they fall short by the +number of such samples. AC, AN and AF are not affected. + +A sample with two different ALT alleles (`GT=1/2`) is counted as heterozygous on +both lines. + ### Warnings afquery emits a `AfqueryWarning` to stderr when a query may silently return fewer or no results. Common causes: @@ -126,7 +139,7 @@ When using `afquery annotate`, the following INFO fields are added to each varia | `AFQUERY_AF` | A (per ALT) | Allele frequency — one value per ALT allele | | `AFQUERY_N_HET` | A (per ALT) | Heterozygous sample count per ALT allele | | `AFQUERY_N_HOM_ALT` | A (per ALT) | Homozygous alt sample count per ALT allele | -| `AFQUERY_N_HOM_REF` | A (per ALT) | Homozygous ref sample count per ALT allele | +| `AFQUERY_N_HOM_REF` | A (per ALT) | Homozygous ref sample count per ALT allele, excluding carriers of other ALT alleles at the position | | `AFQUERY_N_FAIL` | 1 (per site) | Fail sample count — shared across all ALT alleles | | `AFQUERY_N_NO_COVERAGE` | A (per ALT) | Eligible samples whose tech lacks coverage evidence at this position. Always `0` unless a coverage-evidence filter is active. See [Coverage Evidence](../advanced/coverage-evidence.md). | diff --git a/docs/guides/annotate-vcf.md b/docs/guides/annotate-vcf.md index 01b6df5..e9d3b4a 100644 --- a/docs/guides/annotate-vcf.md +++ b/docs/guides/annotate-vcf.md @@ -25,7 +25,7 @@ afquery annotate \ | `AFQUERY_AF` | Float | A (per ALT) | Allele frequency (`AC / AN`) | | `AFQUERY_N_HET` | Integer | A (per ALT) | Heterozygous sample count | | `AFQUERY_N_HOM_ALT` | Integer | A (per ALT) | Homozygous alt sample count | -| `AFQUERY_N_HOM_REF` | Integer | A (per ALT) | Homozygous ref sample count | +| `AFQUERY_N_HOM_REF` | Integer | A (per ALT) | Homozygous ref sample count, excluding carriers of other ALT alleles at the position (including alleles stored in the database but absent from the input record) | | `AFQUERY_N_FAIL` | Integer | 1 (per site) | Eligible samples whose call had FILTER≠PASS. Excluded from AC, but still counted in AN. Mutually exclusive with N_HET/N_HOM_ALT/N_HOM_REF. | | `AFQUERY_N_NO_COVERAGE` | Integer | A (per ALT) | Eligible samples whose tech lacks coverage evidence at this position. Excluded from `N_HOM_REF` to keep AC/AN conservative. Always `0` unless a coverage-evidence filter is active. See [Coverage Evidence](../advanced/coverage-evidence.md). | diff --git a/docs/guides/query.md b/docs/guides/query.md index e0958be..5b7d064 100644 --- a/docs/guides/query.md +++ b/docs/guides/query.md @@ -155,7 +155,7 @@ that fall below a threshold are reported in **N_NO_COVERAGE** instead of N_HOM_R | Flag | Meaning | |------|---------| -| `--min-pass K` | A partially-covered tech is valid for hom-ref at a position only if it has ≥K PASS carriers (het\|hom). Otherwise its non-carrier samples move to `N_NO_COVERAGE`. | +| `--min-pass K` | A partially-covered tech is valid for hom-ref at a position only if it has ≥K PASS carriers (het\|hom) of any ALT allele there. Otherwise its non-carrier samples move to `N_NO_COVERAGE`. | | `--min-observed K` | Same as `--min-pass`, but counts any VCF entry (`het\|hom\|fail`). Useful when you want to include calls that failed FILTER as evidence the position was sequenced. | | `--min-quality-evidence K` | Requires ≥K quality-passing carriers per partially-covered tech. Requires a database built with `--min-dp`, `--min-gq`, `--min-qual`, or `--min-covered`. | @@ -167,8 +167,11 @@ afquery query --db ./db/ --locus chr1:925952 --min-pass 1 afquery query --db ./db/ --region chr1:900000-1000000 --min-observed 2 --min-pass 1 ``` -The genotype invariant becomes: +At a biallelic site the genotype invariant becomes: `N_HET + N_HOM_ALT + N_HOM_REF + N_FAIL + N_NO_COVERAGE = n_eligible`. +At a multi-allelic site, eligible samples that carry only another ALT allele at +the position fall in none of these categories for this allele, so the sum is +lower than `n_eligible` by exactly that number of samples. Fully-covered samples (those whose tech was registered without a BED) are never affected. Carrier samples (het/hom/fail) are never moved to diff --git a/docs/reference/python-api.md b/docs/reference/python-api.md index 3937e12..f4e2840 100644 --- a/docs/reference/python-api.md +++ b/docs/reference/python-api.md @@ -457,8 +457,11 @@ class QueryResult: N_NO_COVERAGE: int # Eligible samples whose tech lacks evidence (excluded from N_HOM_REF) ``` -The new genotype invariant is +At a biallelic site the genotype categories partition the eligible samples: `N_HET + N_HOM_ALT + N_HOM_REF + N_FAIL + N_NO_COVERAGE == n_samples_eligible`. +At a multi-allelic site, eligible samples that carry only another ALT allele at +the position fall in none of these categories for this allele, so the sum is +lower than `n_samples_eligible` by exactly that number of samples. See [Coverage Evidence](../advanced/coverage-evidence.md) for details on `N_NO_COVERAGE`. diff --git a/src/afquery/annotate.py b/src/afquery/annotate.py index 48e501f..4a525ad 100644 --- a/src/afquery/annotate.py +++ b/src/afquery/annotate.py @@ -7,7 +7,6 @@ from .bitmaps import deserialize from .constants import normalize_chrom from .models import AfqueryWarning, SampleFilter -from .ploidy import split_ploidy logger = logging.getLogger(__name__) @@ -33,8 +32,6 @@ def _compute_chunk_annotations( AfqueryWarning, stacklevel=2, ) sample_bm = engine._build_sample_bitmap(sf) - male_bm = engine._male_bm - female_bm = engine._female_bm unique_positions = list({pos for pos, _ref, _alts in records}) @@ -83,6 +80,9 @@ def _compute_chunk_annotations( if pos in valid_pos_set: variant_data[(pos, ref, alt)] = tuple(bytes(b) for b in row[3:3 + n_bitmap_cols]) + unpacked = {key: engine._unpack_bitmaps(raw) for key, raw in variant_data.items()} + sites = engine._site_evidence_by_pos([key[0] for key in unpacked], list(unpacked.values())) + 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] @@ -92,38 +92,19 @@ def _compute_chunk_annotations( continue # dedup if AN == 0: result[key] = (0, 0, False, 0, 0, 0, 0, 0) - elif key in variant_data: - het_bm, hom_bm, fail_bm, filtered_bm, quality_pass_bm = engine._unpack_bitmaps( - variant_data[key] - ) - haploid_elig, diploid_elig = split_ploidy( - eligible, male_bm, female_bm, chrom, pos, engine._genome_build - ) - het_elig = het_bm & eligible - hom_elig = hom_bm & eligible - AC = (len((het_elig | hom_elig) & haploid_elig) - + len(het_elig & diploid_elig) - + 2 * len(hom_elig & 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: int = len(fail_bm & eligible) - no_cov_bm = engine._compute_no_coverage_bm( - eligible, het_bm, hom_bm, fail_bm, - sf.min_pass, sf.min_observed, - filtered_bm=filtered_bm, - quality_pass_bm=quality_pass_bm, - min_quality_evidence=sf.min_quality_evidence, + 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, ) - N_NO_COVERAGE = len(no_cov_bm) - N_HOM_REF = len(eligible) - N_HET - N_HOM_ALT - N_FAIL - N_NO_COVERAGE - result[key] = (AC, AN, True, N_FAIL, N_HET, N_HOM_ALT, N_HOM_REF, N_NO_COVERAGE) + result[key] = (s.AC, AN, True, s.N_FAIL, s.N_HET, s.N_HOM_ALT, + s.N_HOM_REF, s.N_NO_COVERAGE) else: - # Position covered (AN>0) but variant not in Parquet → assume hom-ref - # for all eligible samples. Phase 1/2 filters do not apply because - # there are no carriers at all to evaluate against. - result[key] = (0, AN, False, 0, 0, 0, len(eligible), 0) + # 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. + n_other = len(sites[pos].carrier_bm & eligible) if pos in sites else 0 + result[key] = (0, AN, False, 0, 0, 0, len(eligible) - n_other, 0) return result @@ -170,7 +151,8 @@ def annotate_vcf( }) vcf.add_info_to_header({ "ID": "AFQUERY_N_HOM_REF", "Number": "A", "Type": "Integer", - "Description": "Count of homozygous ref eligible samples per alt allele", + "Description": "Count of homozygous ref eligible samples per alt allele, " + "excluding carriers of other alt alleles at the position", }) vcf.add_info_to_header({ "ID": "AFQUERY_N_FAIL", "Number": "1", "Type": "Integer", diff --git a/src/afquery/dump.py b/src/afquery/dump.py index 51262df..2d062fe 100644 --- a/src/afquery/dump.py +++ b/src/afquery/dump.py @@ -11,7 +11,6 @@ from .bitmaps import deserialize from .constants import normalize_chrom, ALL_CHROMS from .models import SampleFilter -from .ploidy import split_ploidy logger = logging.getLogger(__name__) @@ -180,10 +179,12 @@ def _dump_bucket_worker( pos_cache: dict[int, tuple] = {} group_pos_cache: dict[tuple, tuple] = {} + unpacked = [engine._unpack_bitmaps(row[3:]) for row in rows] + sites = engine._site_evidence_by_pos([row[0] for row in rows], unpacked) + result_rows = [] - for row in rows: + for row, bitmaps in zip(rows, unpacked): pos, ref, alt = row[0], row[1], row[2] - het_bm, hom_bm, fail_bm, filtered_bm, quality_pass_bm = engine._unpack_bitmaps(row[3:]) # Base eligible / AN if pos not in pos_cache: @@ -193,49 +194,27 @@ def _dump_bucket_worker( if AN == 0: continue - haploid_elig, diploid_elig = split_ploidy( - eligible, engine._male_bm, engine._female_bm, chrom, pos, engine._genome_build - ) - het_elig = het_bm & eligible - hom_elig = hom_bm & eligible - AC = ( - len((het_elig | hom_elig) & haploid_elig) - + len(het_elig & diploid_elig) - + 2 * len(hom_elig & diploid_elig) + s = engine._variant_stats( + chrom, pos, eligible, bitmaps, sites[pos], + base_sf.min_pass, base_sf.min_observed, base_sf.min_quality_evidence, ) - if AC == 0 and not include_ac_zero: + if s.AC == 0 and not include_ac_zero: continue # main row filter - N_HET = len(het_elig & diploid_elig) - N_HOM_ALT = ( - len(hom_elig & diploid_elig) + len((het_elig | hom_elig) & haploid_elig) - ) - AF = AC / AN - N_FAIL = len(fail_bm & eligible) - no_cov_bm = engine._compute_no_coverage_bm( - eligible, het_bm, hom_bm, fail_bm, - base_sf.min_pass, base_sf.min_observed, - filtered_bm=filtered_bm, - quality_pass_bm=quality_pass_bm, - min_quality_evidence=base_sf.min_quality_evidence, - ) - N_NO_COVERAGE = len(no_cov_bm) - N_HOM_REF = len(eligible) - N_HET - N_HOM_ALT - N_FAIL - N_NO_COVERAGE - out_row: dict = { "chrom": chrom, "pos": pos, "ref": ref, "alt": alt, - "AC": AC, + "AC": s.AC, "AN": AN, - "AF": AF, - "N_HET": N_HET, - "N_HOM_ALT": N_HOM_ALT, - "N_HOM_REF": N_HOM_REF, - "N_FAIL": N_FAIL, - "N_NO_COVERAGE": N_NO_COVERAGE, + "AF": s.AC / AN, + "N_HET": s.N_HET, + "N_HOM_ALT": s.N_HOM_ALT, + "N_HOM_REF": s.N_HOM_REF, + "N_FAIL": s.N_FAIL, + "N_NO_COVERAGE": s.N_NO_COVERAGE, } # Per-group columns @@ -246,42 +225,19 @@ def _dump_bucket_worker( group_pos_cache[cache_key] = engine._compute_eligible(chrom, pos, g_bm) g_eligible, g_AN = group_pos_cache[cache_key] - g_haploid, g_diploid = split_ploidy( - g_eligible, engine._male_bm, engine._female_bm, chrom, pos, engine._genome_build - ) - g_het_elig = het_bm & g_eligible - g_hom_elig = hom_bm & g_eligible - g_AC = ( - len((g_het_elig | g_hom_elig) & g_haploid) - + len(g_het_elig & g_diploid) - + 2 * len(g_hom_elig & g_diploid) - ) - g_N_HET = len(g_het_elig & g_diploid) - g_N_HOM_ALT = ( - len(g_hom_elig & g_diploid) + len((g_het_elig | g_hom_elig) & g_haploid) - ) - g_AF = g_AC / g_AN if g_AN > 0 else 0.0 - g_N_FAIL = len(fail_bm & g_eligible) - g_no_cov_bm = engine._compute_no_coverage_bm( - g_eligible, het_bm, hom_bm, fail_bm, - g_sf.min_pass, g_sf.min_observed, - filtered_bm=filtered_bm, - quality_pass_bm=quality_pass_bm, - min_quality_evidence=g_sf.min_quality_evidence, - ) - g_N_NO_COVERAGE = len(g_no_cov_bm) - g_N_HOM_REF = ( - len(g_eligible) - g_N_HET - g_N_HOM_ALT - g_N_FAIL - g_N_NO_COVERAGE + g = engine._variant_stats( + chrom, pos, g_eligible, bitmaps, sites[pos], + g_sf.min_pass, g_sf.min_observed, g_sf.min_quality_evidence, ) - out_row[f"AC_{label}"] = g_AC + out_row[f"AC_{label}"] = g.AC out_row[f"AN_{label}"] = g_AN - out_row[f"AF_{label}"] = g_AF - 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 - out_row[f"N_FAIL_{label}"] = g_N_FAIL - out_row[f"N_NO_COVERAGE_{label}"] = g_N_NO_COVERAGE + out_row[f"AF_{label}"] = g.AC / g_AN if g_AN > 0 else 0.0 + 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 + out_row[f"N_FAIL_{label}"] = g.N_FAIL + out_row[f"N_NO_COVERAGE_{label}"] = g.N_NO_COVERAGE result_rows.append(out_row) diff --git a/src/afquery/query.py b/src/afquery/query.py index b2aac57..c383963 100644 --- a/src/afquery/query.py +++ b/src/afquery/query.py @@ -2,6 +2,7 @@ import sqlite3 import warnings from pathlib import Path +from typing import NamedTuple import duckdb from pyroaring import BitMap @@ -16,6 +17,22 @@ BATCH_IN_THRESHOLD = 10_000 +class SiteEvidence(NamedTuple): + """Carriers of any allele at one position, pooled across its Parquet rows.""" + pass_bm: BitMap + carrier_bm: BitMap + quality_pass_bm: "BitMap | None" + + +class VariantStats(NamedTuple): + AC: int + N_HET: int + N_HOM_ALT: int + N_HOM_REF: int + N_FAIL: int + N_NO_COVERAGE: int + + def _parse_schema_version(s: str) -> tuple[int, ...]: try: return tuple(int(p) for p in s.split(".")) @@ -203,12 +220,14 @@ def _compute_no_coverage_bm( filtered_bm: "BitMap | None" = None, quality_pass_bm: "BitMap | None" = None, min_quality_evidence: int = 0, + site: "SiteEvidence | None" = None, ) -> BitMap: """Return bitmap of WES non-carrier samples that should not be assumed hom-ref. Phase 1 (query-time): count-based per-tech gate using existing bitmaps. Phase 2 (build-time): stored filtered_bitmap + optional quality_pass gate. - Results are unioned; carriers (het/hom/fail) are never included. + Results are unioned; carriers (het/hom/fail) are never included, nor are + carriers of any other allele at the position when ``site`` is given. """ if min_quality_evidence > 0 and not self._has_coverage_data: raise ValueError( @@ -216,7 +235,15 @@ def _compute_no_coverage_bm( "Re-create with --min-dp / --min-gq to use --min-quality-evidence." ) no_cov = BitMap() + pass_carriers = het_bm | hom_bm all_carriers = het_bm | hom_bm | fail_bm + # Coverage evidence is a property of the position: a call for any + # allele there shows it was sequenced, and that sample is not uncovered. + if site is not None: + pass_carriers = pass_carriers | site.pass_bm + all_carriers = all_carriers | site.carrier_bm + if site.quality_pass_bm is not None: + quality_pass_bm = site.quality_pass_bm # Phase 1: count-based gate if min_pass > 0 or min_observed > 0: @@ -225,14 +252,15 @@ def _compute_no_coverage_bm( tech_eligible = eligible & tech_bm if len(tech_eligible) == 0: continue - pass_count = len((het_bm | hom_bm) & tech_eligible) + pass_count = len(pass_carriers & tech_eligible) observed_count = len(all_carriers & tech_eligible) if pass_count < min_pass or observed_count < min_observed: no_cov |= tech_eligible - (all_carriers & tech_eligible) - # Phase 2: stored filtered_bitmap + # Phase 2: stored filtered_bitmap (built per allele, so it may hold + # carriers of another allele at the same position) if filtered_bm is not None: - no_cov |= filtered_bm & eligible + no_cov |= (filtered_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: @@ -282,6 +310,86 @@ def _unpack_bitmaps( quality_pass_bm = None return het_bm, hom_bm, fail_bm, filtered_bm, quality_pass_bm + @staticmethod + 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() + 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 + 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) + + @classmethod + def _site_evidence_by_pos( + cls, positions: list[int], unpacked_rows: list[tuple], + ) -> dict[int, SiteEvidence]: + """Group unpacked rows by position and pool each group.""" + by_pos: dict[int, list[tuple]] = {} + for pos, bitmaps in zip(positions, unpacked_rows): + by_pos.setdefault(pos, []).append(bitmaps) + return {pos: cls._site_evidence(group) for pos, group in by_pos.items()} + + def _variant_stats( + self, + chrom: str, + pos: int, + eligible: BitMap, + bitmaps: tuple, + site: SiteEvidence, + min_pass: int, + min_observed: int, + min_quality_evidence: int, + ) -> VariantStats: + """Genotype tallies for one allele among the eligible samples. + + Eligible samples carrying only another allele at this position 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 + haploid_elig, diploid_elig = split_ploidy( + eligible, self._male_bm, self._female_bm, chrom, pos, self._genome_build + ) + het_elig = het_bm & eligible + hom_elig = hom_bm & eligible + AC = (len((het_elig | hom_elig) & haploid_elig) + + len(het_elig & diploid_elig) + + 2 * len(hom_elig & 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) + no_cov_bm = self._compute_no_coverage_bm( + eligible, het_bm, hom_bm, fail_bm, + min_pass, min_observed, + filtered_bm=filtered_bm, + quality_pass_bm=quality_pass_bm, + min_quality_evidence=min_quality_evidence, + site=site, + ) + N_NO_COVERAGE = len(no_cov_bm) + other_allele = (site.carrier_bm - (het_bm | hom_bm | fail_bm)) & eligible + N_HOM_REF = (len(eligible) - N_HET - N_HOM_ALT - N_FAIL - N_NO_COVERAGE + - len(other_allele)) + return VariantStats(AC, 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, + ) -> 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, + 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, + ) + def _compute_eligible( self, chrom: str, @@ -364,38 +472,16 @@ def query(self, params: QueryParams) -> list[QueryResult]: if not rows: return [] - results = [] sf = params.filter - for row in rows: - ref, alt = row[0], row[1] - het_bm, hom_bm, fail_bm, filtered_bm, quality_pass_bm = self._unpack_bitmaps(row[2:]) - haploid_elig, diploid_elig = split_ploidy( - eligible, self._male_bm, self._female_bm, chrom, pos, self._genome_build - ) - het_elig = het_bm & eligible - hom_elig = hom_bm & eligible - AC = (len((het_elig | hom_elig) & haploid_elig) - + len(het_elig & diploid_elig) - + 2 * len(hom_elig & diploid_elig)) - N_HET = len(het_elig & diploid_elig) - N_HOM_ALT = len(hom_elig & diploid_elig) + len((het_elig | hom_elig) & haploid_elig) - AF = AC / AN if AN > 0 else None - N_FAIL = len(fail_bm & eligible) - no_cov_bm = self._compute_no_coverage_bm( - eligible, het_bm, hom_bm, fail_bm, - sf.min_pass, sf.min_observed, - filtered_bm=filtered_bm, - quality_pass_bm=quality_pass_bm, - min_quality_evidence=sf.min_quality_evidence, + unpacked = [(row[0], row[1], self._unpack_bitmaps(row[2:])) for row in rows] + site = self._site_evidence([bitmaps for _ref, _alt, bitmaps in unpacked]) + results = [] + for ref, alt, bitmaps in unpacked: + stats = self._variant_stats( + chrom, pos, eligible, bitmaps, site, + sf.min_pass, sf.min_observed, sf.min_quality_evidence, ) - N_NO_COVERAGE = len(no_cov_bm) - N_HOM_REF = len(eligible) - N_HET - N_HOM_ALT - N_FAIL - N_NO_COVERAGE - results.append(QueryResult( - variant=VariantKey(chrom=chrom, pos=pos, ref=ref, alt=alt), - AC=AC, AN=AN, AF=AF, n_samples_eligible=len(eligible), - N_HET=N_HET, N_HOM_ALT=N_HOM_ALT, N_HOM_REF=N_HOM_REF, - N_FAIL=N_FAIL, N_NO_COVERAGE=N_NO_COVERAGE, - )) + results.append(self._make_result(chrom, pos, ref, alt, AN, 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: @@ -566,40 +652,18 @@ def _query_region_inner( eligible, AN = self._compute_eligible(chrom, pos, sample_bm) pos_data[pos] = (eligible, AN) + rows = [row for row in rows if pos_data[row[0]][1] > 0] + unpacked = [self._unpack_bitmaps(row[3:]) for row in rows] + sites = self._site_evidence_by_pos([row[0] for row in rows], unpacked) results = [] - for row in rows: + for row, bitmaps in zip(rows, unpacked): pos, ref, alt = row[0], row[1], row[2] eligible, AN = pos_data[pos] - if AN == 0: - continue - het_bm, hom_bm, fail_bm, filtered_bm, quality_pass_bm = self._unpack_bitmaps(row[3:]) - haploid_elig, diploid_elig = split_ploidy( - eligible, self._male_bm, self._female_bm, chrom, pos, self._genome_build + stats = self._variant_stats( + chrom, pos, eligible, bitmaps, sites[pos], + min_pass, min_observed, min_quality_evidence, ) - het_elig = het_bm & eligible - hom_elig = hom_bm & eligible - AC = (len((het_elig | hom_elig) & haploid_elig) - + len(het_elig & diploid_elig) - + 2 * len(hom_elig & diploid_elig)) - N_HET = len(het_elig & diploid_elig) - N_HOM_ALT = len(hom_elig & diploid_elig) + len((het_elig | hom_elig) & haploid_elig) - AF = AC / AN if AN > 0 else None - N_FAIL = len(fail_bm & eligible) - no_cov_bm = self._compute_no_coverage_bm( - eligible, het_bm, hom_bm, fail_bm, - min_pass, min_observed, - filtered_bm=filtered_bm, - quality_pass_bm=quality_pass_bm, - min_quality_evidence=min_quality_evidence, - ) - N_NO_COVERAGE = len(no_cov_bm) - N_HOM_REF = len(eligible) - N_HET - N_HOM_ALT - N_FAIL - N_NO_COVERAGE - results.append(QueryResult( - variant=VariantKey(chrom=chrom, pos=pos, ref=ref, alt=alt), - AC=AC, AN=AN, AF=AF, n_samples_eligible=len(eligible), - N_HET=N_HET, N_HOM_ALT=N_HOM_ALT, N_HOM_REF=N_HOM_REF, - N_FAIL=N_FAIL, N_NO_COVERAGE=N_NO_COVERAGE, - )) + results.append(self._make_result(chrom, pos, ref, alt, AN, eligible, stats)) results.sort(key=lambda r: (r.variant.pos, r.variant.alt)) return results @@ -658,40 +722,21 @@ def _query_batch_inner( ).fetchall() con.close() + # Pooled before the request filter: an allele nobody asked for still + # has carriers who must not be counted as hom-ref for the others. + unpacked = [self._unpack_bitmaps(row[3:]) for row in rows] + sites = self._site_evidence_by_pos([row[0] for row in rows], unpacked) results = [] - for row in rows: + for row, bitmaps in zip(rows, unpacked): pos, ref, alt = row[0], row[1], row[2] if (pos, ref, alt) not in requested_variants: continue eligible, AN = pos_data[pos] - het_bm, hom_bm, fail_bm, filtered_bm, quality_pass_bm = self._unpack_bitmaps(row[3:]) - haploid_elig, diploid_elig = split_ploidy( - eligible, self._male_bm, self._female_bm, chrom, pos, self._genome_build - ) - het_elig = het_bm & eligible - hom_elig = hom_bm & eligible - AC = (len((het_elig | hom_elig) & haploid_elig) - + len(het_elig & diploid_elig) - + 2 * len(hom_elig & diploid_elig)) - N_HET = len(het_elig & diploid_elig) - N_HOM_ALT = len(hom_elig & diploid_elig) + len((het_elig | hom_elig) & haploid_elig) - AF = AC / AN if AN > 0 else None - N_FAIL = len(fail_bm & eligible) - no_cov_bm = self._compute_no_coverage_bm( - eligible, het_bm, hom_bm, fail_bm, - min_pass, min_observed, - filtered_bm=filtered_bm, - quality_pass_bm=quality_pass_bm, - min_quality_evidence=min_quality_evidence, + stats = self._variant_stats( + chrom, pos, eligible, bitmaps, sites[pos], + min_pass, min_observed, min_quality_evidence, ) - N_NO_COVERAGE = len(no_cov_bm) - N_HOM_REF = len(eligible) - N_HET - N_HOM_ALT - N_FAIL - N_NO_COVERAGE - results.append(QueryResult( - variant=VariantKey(chrom=chrom, pos=pos, ref=ref, alt=alt), - AC=AC, AN=AN, AF=AF, n_samples_eligible=len(eligible), - N_HET=N_HET, N_HOM_ALT=N_HOM_ALT, N_HOM_REF=N_HOM_REF, - N_FAIL=N_FAIL, N_NO_COVERAGE=N_NO_COVERAGE, - )) + results.append(self._make_result(chrom, pos, ref, alt, AN, eligible, stats)) results.sort(key=lambda r: (r.variant.pos, r.variant.alt)) return results @@ -764,6 +809,7 @@ def variant_info(self, params: QueryParams) -> list[SampleCarrier]: if not rows: return [] + site = self._site_evidence([self._unpack_bitmaps(r[3:]) for r in rows]) if params.ref is not None: rows = [r for r in rows if r[1] == params.ref] if params.alt is not None: @@ -798,6 +844,7 @@ def variant_info(self, params: QueryParams) -> list[SampleCarrier]: filtered_bm=filtered_bm, quality_pass_bm=quality_pass_bm, min_quality_evidence=sf.min_quality_evidence, + site=site, ) no_cov_elig = no_cov_bm & eligible diff --git a/tests/oracle.py b/tests/oracle.py index 5cfc614..e77d740 100644 --- a/tests/oracle.py +++ b/tests/oracle.py @@ -162,13 +162,25 @@ def expect(self, chrom, pos, ref, alt, samples=None) -> dict: AC += 2 n_hom_alt += 1 + # A sample carrying only another allele at this position is not hom-ref + # for this one: at a multi-allelic site the tallies do not add up to + # n_eligible. + n_other = sum( + 1 for s in elig + if 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) + ) + ) + return { "AC": AC, "AN": AN, "N_HET": n_het, "N_HOM_ALT": n_hom_alt, "N_FAIL": n_fail, - "N_HOM_REF": len(elig) - n_het - n_hom_alt - n_fail, + "N_HOM_REF": len(elig) - n_het - n_hom_alt - n_fail - n_other, } def variants(self) -> list[tuple[str, int, str, str]]: diff --git a/tests/test_multiallelic.py b/tests/test_multiallelic.py new file mode 100644 index 0000000..79344b6 --- /dev/null +++ b/tests/test_multiallelic.py @@ -0,0 +1,153 @@ +"""Genotype tallies at multi-allelic sites. + +Each ALT allele at a position is its own row. A sample carrying only another +allele there is neither a carrier nor hom-ref for this one, so it is left out of +every tally and N_HET + N_HOM_ALT + N_HOM_REF + N_FAIL + N_NO_COVERAGE falls +short of n_eligible by exactly that number of samples. +""" +import csv +import io +import warnings + +import cyvcf2 +import pytest + +from afquery.models import AfqueryWarning + +from test_variant_info import _build_db + +SITE = ("chr1", 5000) + +COHORT = [ + ("A_HET", "female", "WGS", [("chr1", 5000, "G", "A", "0/1")]), + ("A_HOM", "male", "WGS", [("chr1", 5000, "G", "A", "1/1")]), + ("T_HET", "female", "WGS", [("chr1", 5000, "G", "T", "0/1")]), + ("AT", "male", "WGS", [("chr1", 5000, "G", "A,T", "1/2")]), + ("REF1", "female", "WGS", [("chr1", 100, "C", "G", "0/1")]), + ("REF2", "male", "WGS", [("chr1", 100, "C", "G", "0/1")]), +] + +# (N_HET, N_HOM_ALT, N_HOM_REF, samples carrying only the other allele) +EXPECTED = { + "A": (2, 1, 2, 1), # het A_HET, AT; hom A_HOM; hom-ref REF1, REF2; other T_HET + "T": (2, 0, 2, 2), # het T_HET, AT; hom-ref REF1, REF2; other A_HET, A_HOM +} + + +@pytest.fixture(scope="module") +def db(tmp_path_factory): + return _build_db(tmp_path_factory.mktemp("multiallelic"), COHORT) + + +def _tallies(r): + return (r.N_HET, r.N_HOM_ALT, r.N_HOM_REF) + + +def _check(results): + by_alt = {r.variant.alt: r for r in results if (r.variant.chrom, r.variant.pos) == SITE} + assert set(by_alt) == set(EXPECTED) + for alt, (het, hom, hom_ref, other) in EXPECTED.items(): + r = by_alt[alt] + assert _tallies(r) == (het, hom, hom_ref), alt + assert r.n_samples_eligible == 6 + total = r.N_HET + r.N_HOM_ALT + r.N_HOM_REF + r.N_FAIL + r.N_NO_COVERAGE + assert total == r.n_samples_eligible - other, alt + + +def test_point_query(db): + _check(db.query(*SITE)) + + +def test_region_query(db): + _check(db.query_region("chr1", 1, 10000)) + + +def test_batch_query_one_allele_requested(db): + """The allele not asked for still keeps its carriers out of hom-ref.""" + [r] = db.query_batch("chr1", [(5000, "G", "A")]) + assert _tallies(r) == EXPECTED["A"][:3] + + +def test_batch_multi_query(db): + results = db.query_batch_multi([("chr1", 5000, "G", "T"), ("chr1", 5000, "G", "A")]) + _check(results) + + +def test_dump_matches_query(db): + buf = io.StringIO() + db.dump(output=buf, by_sex=True) + rows =[r for r in csv.DictReader(io.StringIO(buf.getvalue())) if r["pos"] == "5000"] + assert {r["alt"] for r in rows} == set(EXPECTED) + for row in rows: + [base] = db.query(*SITE, alt=row["alt"]) + assert int(row["N_HOM_REF"]) == base.N_HOM_REF == EXPECTED[row["alt"]][2] + for sex in ("male", "female"): + [g] = db.query(*SITE, alt=row["alt"], sex=sex) + assert int(row[f"N_HOM_REF_{sex}"]) == g.N_HOM_REF + assert int(row[f"N_HET_{sex}"]) == g.N_HET + + +def test_annotate_matches_query(db, tmp_path): + """C is not stored at the site: its hom-ref excludes carriers of 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" + db.annotate_vcf(str(vcf_in), str(out)) + [v] = list(cyvcf2.VCF(str(out))) + assert tuple(v.INFO.get("AFQUERY_N_HOM_REF")) == (2, 2, 2) + assert tuple(v.INFO.get("AFQUERY_N_HET")) == (2, 2, 0) + + +def test_variant_info_agrees(db): + with warnings.catch_warnings(): + warnings.simplefilter("ignore", AfqueryWarning) + carriers = db.variant_info(*SITE) + names = sorted(c.sample_name for c in carriers) + assert names == ["AT", "AT", "A_HET", "A_HOM", "T_HET"] + + +# --------------------------------------------------------------------------- +# Coverage evidence is judged per position, not per allele +# --------------------------------------------------------------------------- + +PANEL_COHORT = COHORT + [ + ("P_T", "female", "PANEL", [("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")]), +] + + +@pytest.fixture(scope="module") +def panel_db(tmp_path_factory): + return _build_db( + tmp_path_factory.mktemp("multiallelic_panel"), PANEL_COHORT, + beds={"PANEL": "chr1\t0\t6000\n"}, + ) + + +def test_min_pass_counts_any_allele_as_evidence(panel_db): + """The panel has no A carrier but one T carrier: it was sequenced here, so + with --min-pass 1 its non-carriers stay hom-ref for A as well.""" + by_alt = {r.variant.alt: r for r in panel_db.query(*SITE, min_pass=1)} + assert by_alt["A"].N_NO_COVERAGE == 0 + assert by_alt["A"].N_HOM_REF == 4 # REF1, REF2, P_R1, P_R2 + assert by_alt["T"].N_NO_COVERAGE == 0 + + +def test_min_pass_failure_never_marks_carriers_uncovered(panel_db): + """With --min-pass 2 the panel fails the gate; P_T has a call, so only the + panel's true non-carriers move to N_NO_COVERAGE, for every allele.""" + by_alt = {r.variant.alt: r for r in panel_db.query(*SITE, min_pass=2)} + for alt in ("A", "T"): + assert by_alt[alt].N_NO_COVERAGE == 2, alt # P_R1, P_R2 + assert by_alt["A"].N_HOM_REF == 2 # REF1, REF2 + with warnings.catch_warnings(): + warnings.simplefilter("ignore", AfqueryWarning) + 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"] diff --git a/tests/test_oracle_consistency.py b/tests/test_oracle_consistency.py index 121237a..d7b97ab 100644 --- a/tests/test_oracle_consistency.py +++ b/tests/test_oracle_consistency.py @@ -162,3 +162,62 @@ def test_after_add_samples_new_phenotype_matches_oracle(grown_db, data_dir, tmp_ cohort = oracle.Cohort(grow_manifest, bed_dir=data_dir / "beds") _assert_matches(db, cohort, samples=["S10", "S11"], phenotype=["ZZNEW"], label="added-only") + + +# --------------------------------------------------------------------------- +# Multi-allelic cohort +# --------------------------------------------------------------------------- + +def test_multiallelic_cohort_matches_oracle(tmp_path): + """Two ALT alleles at one position, a 1/2 carrier, partial capture and a + FILTER failure: every allele's tallies agree with the oracle.""" + calls = { + # name: (sex, tech, [(chrom, pos, ref, alt, gt, filter)]) + "A_HET": ("female", "WGS", [("chr1", 5000, "G", "A", "0/1", "PASS")]), + "A_HOM": ("male", "WGS", [("chr1", 5000, "G", "A", "1/1", "PASS")]), + "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")]), + "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")]), + "P_OFF": ("female", "PANEL", [("chr1", 9000, "C", "G", "0/1", "PASS"), + ("chr1", 9000, "C", "T", "0/1", "PASS")]), + "X_A": ("male", "WGS", [("chrX", 5000000, "A", "G,C", "1", "PASS")]), + "X_C": ("female", "WGS", [("chrX", 5000000, "A", "C", "0/1", "PASS")]), + } + beds = tmp_path / "beds" + beds.mkdir() + (beds / "PANEL.bed").write_text("chr1\t0\t6000\n") + manifest = ["sample_name\tsex\ttech_name\tvcf_path\tphenotype_codes"] + for name, (sex, tech, records) in calls.items(): + vcf = tmp_path / f"{name}.vcf" + with open(vcf, "w") as f: + f.write("##fileformat=VCFv4.2\n") + f.write('##FILTER=\n') + for contig in sorted({r[0] for r in records}): + f.write(f"##contig=\n") + f.write('##FORMAT=\n') + f.write(f"#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t{name}\n") + for chrom, pos, ref, alt, gt, flt in records: + f.write(f"{chrom}\t{pos}\t.\t{ref}\t{alt}\t.\t{flt}\t.\tGT\t{gt}\n") + manifest.append(f"{name}\t{sex}\t{tech}\t{vcf}\tE11.9") + manifest_path = tmp_path / "manifest.tsv" + manifest_path.write_text("\n".join(manifest) + "\n") + + db = tmp_path / "db" + run_preprocess( + manifest_path=str(manifest_path), output_dir=str(db), + genome_build="GRCh37", bed_dir=str(beds), threads=1, + ) + 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 + # 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) From 672640319cc84e1879a4c624e95546987411b0c9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20L=C3=B3pez=20L=C3=B3pez?= Date: Wed, 30 Sep 2026 12:59:53 +0200 Subject: [PATCH 3/5] perf(query): skip allele pooling at single-allele positions Most positions hold one allele, where the pooled evidence equals the row's own bitmaps. Pass None instead of building it, and index stored alleles by position in annotate rather than scanning them per record. --- src/afquery/annotate.py | 12 ++++++++++-- src/afquery/query.py | 32 +++++++++++++++++++++----------- 2 files changed, 31 insertions(+), 13 deletions(-) diff --git a/src/afquery/annotate.py b/src/afquery/annotate.py index 4a525ad..ccede51 100644 --- a/src/afquery/annotate.py +++ b/src/afquery/annotate.py @@ -81,7 +81,13 @@ def _compute_chunk_annotations( variant_data[(pos, ref, alt)] = tuple(bytes(b) for b in row[3:3 + n_bitmap_cols]) unpacked = {key: engine._unpack_bitmaps(raw) for key, raw in variant_data.items()} - sites = engine._site_evidence_by_pos([key[0] for key in unpacked], list(unpacked.values())) + stored_by_pos: dict[int, list[tuple]] = {} + for (pos, _ref, _alt), bitmaps in unpacked.items(): + stored_by_pos.setdefault(pos, []).append(bitmaps) + sites = { + pos: engine._site_evidence(group) if len(group) > 1 else None + for pos, group in stored_by_pos.items() + } result: dict[tuple[int, str, str], tuple[int, int, bool, int, int, int, int, int]] = {} for pos, ref, alts in records: @@ -103,7 +109,9 @@ def _compute_chunk_annotations( # 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. - n_other = len(sites[pos].carrier_bm & eligible) if pos in sites else 0 + 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) return result diff --git a/src/afquery/query.py b/src/afquery/query.py index c383963..c8ae6df 100644 --- a/src/afquery/query.py +++ b/src/afquery/query.py @@ -326,12 +326,18 @@ def _site_evidence(unpacked_rows: list[tuple]) -> SiteEvidence: @classmethod def _site_evidence_by_pos( cls, positions: list[int], unpacked_rows: list[tuple], - ) -> dict[int, SiteEvidence]: - """Group unpacked rows by position and pool each group.""" + ) -> dict[int, "SiteEvidence | None"]: + """Group unpacked rows by position and pool each group. + + Single-allele positions map to None: there is nothing to pool. + """ by_pos: dict[int, list[tuple]] = {} for pos, bitmaps in zip(positions, unpacked_rows): by_pos.setdefault(pos, []).append(bitmaps) - return {pos: cls._site_evidence(group) for pos, group in by_pos.items()} + return { + pos: cls._site_evidence(group) if len(group) > 1 else None + for pos, group in by_pos.items() + } def _variant_stats( self, @@ -339,15 +345,16 @@ def _variant_stats( pos: int, eligible: BitMap, bitmaps: tuple, - site: SiteEvidence, + site: "SiteEvidence | None", min_pass: int, min_observed: int, min_quality_evidence: int, ) -> VariantStats: """Genotype tallies for one allele among the eligible samples. - Eligible samples carrying only another allele at this position are in - none of the returned categories, so at a multi-allelic site + ``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 + 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. """ @@ -372,9 +379,10 @@ def _variant_stats( site=site, ) N_NO_COVERAGE = len(no_cov_bm) - other_allele = (site.carrier_bm - (het_bm | hom_bm | fail_bm)) & eligible - N_HOM_REF = (len(eligible) - N_HET - N_HOM_ALT - N_FAIL - N_NO_COVERAGE - - len(other_allele)) + N_OTHER = 0 + if site is not None: + N_OTHER = len((site.carrier_bm - (het_bm | hom_bm | 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) @staticmethod @@ -474,7 +482,8 @@ def query(self, params: QueryParams) -> list[QueryResult]: sf = params.filter unpacked = [(row[0], row[1], self._unpack_bitmaps(row[2:])) for row in rows] - site = self._site_evidence([bitmaps for _ref, _alt, bitmaps in unpacked]) + site = (self._site_evidence([bitmaps for _ref, _alt, bitmaps in unpacked]) + if len(unpacked) > 1 else None) results = [] for ref, alt, bitmaps in unpacked: stats = self._variant_stats( @@ -809,7 +818,8 @@ def variant_info(self, params: QueryParams) -> list[SampleCarrier]: if not rows: return [] - site = self._site_evidence([self._unpack_bitmaps(r[3:]) for r in rows]) + site = (self._site_evidence([self._unpack_bitmaps(r[3:]) for r in rows]) + if len(rows) > 1 else None) if params.ref is not None: rows = [r for r in rows if r[1] == params.ref] if params.alt is not None: From 56476888d5ce6dd068fd25aba64a37b045ee161d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20L=C3=B3pez=20L=C3=B3pez?= Date: Wed, 30 Sep 2026 13:03:35 +0200 Subject: [PATCH 4/5] fix(normalize): split multi-allelic records and stop dropping them as duplicates normalize_vcf.sh ran bcftools norm without -m, so multi-allelic records were kept whole and their indels were not trimmed per allele: the 2:136592357 CAAAAAAA>C,CAAAAAAAAAAAA record stored its insertion as CAAAAAAA>CAAAAAAAAAAAA instead of C>CAAAAA, splitting one variant into two representations across samples. It also passed -d twice; the second value, "both", treats any two indels at one position as duplicates, so once the record is split one of the two alleles was discarded. Use -m -both -d exact. bcftools refuses to combine -m and -d before 1.20, so the script now checks the version and exits with a clear message on older releases. Closes #39 --- docs/getting-started/preprocessing.md | 12 ++++++++++-- resources/normalize_vcf.sh | 13 ++++++++++++- 2 files changed, 22 insertions(+), 3 deletions(-) diff --git a/docs/getting-started/preprocessing.md b/docs/getting-started/preprocessing.md index 43a2f14..08e6f64 100644 --- a/docs/getting-started/preprocessing.md +++ b/docs/getting-started/preprocessing.md @@ -2,7 +2,15 @@ AFQuery ingests single-sample, normalized VCF files. While AFQuery itself does not perform VCF normalization, the accuracy of your allele frequency estimates depends on the quality and consistency of your input VCFs. This page explains why normalization matters and provides a reference pipeline. -A ready-to-use normalization script is provided at [`resources/normalize_vcf.sh`](https://github.com/babelomics/afquery/blob/master/resources/normalize_vcf.sh). +A ready-to-use normalization script is provided at [`resources/normalize_vcf.sh`](https://github.com/babelomics/afquery/blob/master/resources/normalize_vcf.sh). It requires **bcftools 1.20 or later**: older releases refuse to split multi-allelic records and remove duplicates in the same `bcftools norm` call, and the script exits with an error if it finds one. + +!!! warning "VCFs normalized with the script from afquery 0.4.2 or earlier" + Those versions of the script did not split multi-allelic records, and removed + one of two different indels at the same position as a duplicate. Indels from + multi-allelic sites could therefore be stored untrimmed (e.g. `CAAAAAAA>CAAAAAAAAAAAA` + instead of `C>CAAAAA`), so the same variant appeared under two representations, or + were dropped. Re-normalize those VCFs with the current script and rebuild the + database to correct them. --- @@ -17,7 +25,7 @@ Multi-allelic variants and complex indels can be represented in multiple equival - Multi-allelic sites may not be decomposed into biallelic records - Duplicate records can inflate AC -`bcftools norm` left-aligns indels against the reference genome and decomposes multi-allelic sites, ensuring consistent representation across samples. +`bcftools norm -m -both` decomposes multi-allelic sites into biallelic records, and left-aligns and trims each allele against the reference genome, ensuring consistent representation across samples. `-d exact` then drops only records whose alleles are identical, so two different alleles at the same position are both kept. ### Ploidy correction for sex chromosomes diff --git a/resources/normalize_vcf.sh b/resources/normalize_vcf.sh index f3a046c..9b9a5df 100755 --- a/resources/normalize_vcf.sh +++ b/resources/normalize_vcf.sh @@ -25,6 +25,17 @@ REF="$3" GENDER="$4" THREADS="${5:-4}" +# ----------------------------- +# bcftools version (< 1.20 refuses to combine -m and -d in bcftools norm) +# ----------------------------- +MIN_BCFTOOLS="1.20" +BCFTOOLS_VERSION=$(bcftools --version | head -n1 | awk '{print $2}') +if [[ "$(printf '%s\n' "${MIN_BCFTOOLS}" "${BCFTOOLS_VERSION}" | sort -V | head -n1)" != "${MIN_BCFTOOLS}" ]] +then + echo "ERROR: bcftools >= ${MIN_BCFTOOLS} is required (found ${BCFTOOLS_VERSION})" + exit 1 +fi + OUT="${VCF_ID}_norm.vcf.gz" CHR_MAP="chr_list.txt" GENDER_FILE="${VCF_ID}_gender.txt" @@ -88,8 +99,8 @@ bcftools annotate \ --rename-chrs "${CHR_MAP}" \ "${VCF}" | \ bcftools norm \ + -m -both \ -d exact \ - -d both \ -f "${REF}" \ --check-ref ws \ --targets "${TARGETS}" | \ From 5629bf6f682df46b9dd0fcbe85d427ce33374cec Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20L=C3=B3pez=20L=C3=B3pez?= Date: Wed, 30 Sep 2026 13:12:26 +0200 Subject: [PATCH 5/5] docs: match debugging advice to query output and scope --min-covered query prints nothing for a technology that does not cover the position, so the capture-index check now looks for the 'No variants found' message. State that --min-covered is still evaluated per allele at build time, while --min-quality-evidence counts carriers of any allele. --- docs/advanced/coverage-evidence.md | 4 ++-- docs/advanced/debugging-results.md | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/advanced/coverage-evidence.md b/docs/advanced/coverage-evidence.md index 762bc54..9e94460 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`. | +| `--min-covered K`| Per partially-covered tech, the position is "trusted" only if at least K of its carriers pass the quality thresholds. Non-carriers of failing positions are recorded as `N_NO_COVERAGE`. The gate is evaluated at build time for each ALT allele separately; samples carrying another allele at the position are never counted as `N_NO_COVERAGE`. | A carrier counts as quality-passing only if **all** active thresholds hold (unset thresholds are simply ignored). At least one of these flags must be @@ -137,7 +137,7 @@ afquery query --db ./db/ --locus chr1:925952 --min-quality-evidence 5 ``` `--min-quality-evidence K` requires each partially-covered tech to have ≥K -quality-passing carriers at the position. Non-carriers of failing techs +quality-passing carriers of any ALT allele at the position. Non-carriers of failing techs (other than those already filtered at build time) move to `N_NO_COVERAGE`. Running the flag against a database that was not built with quality data diff --git a/docs/advanced/debugging-results.md b/docs/advanced/debugging-results.md index 3de8794..03834c4 100644 --- a/docs/advanced/debugging-results.md +++ b/docs/advanced/debugging-results.md @@ -16,7 +16,7 @@ AN=0 means no eligible samples at the queried position. Work through these check | Position exists in database | `afquery query --db ./db/ --locus chr1:12345678` | If no result at all, the variant was not observed in any sample during ingestion. | | BED coverage (WES) | `afquery info --db ./db/` | If all eligible samples are WES and the position is outside capture regions, AN=0 is correct. | | Sample filter too restrictive | Remove `--phenotype` and `--sex` filters | Query with no filters first. If AN>0 without filters, the filter is excluding all samples. | -| WES samples missing from AN | `afquery query --db ./db/ --locus chr1:12345678 --tech WES_kit_A` | If `n_eligible=0` for a technology whose kit should cover the position, that technology's capture index is not matching it. `afquery` warns on open when a capture BED is empty or matches no known chromosome, naming the offending contigs. | +| WES samples missing from AN | `afquery query --db ./db/ --locus chr1:12345678 --tech WES_kit_A` | If this prints `No variants found for the given filters.` for a technology whose kit should cover the position, while the unfiltered query returns the variant, that technology's capture index is not matching the position. `afquery` warns on open when a capture BED is empty or matches no known chromosome, naming the offending contigs. | | AN=0 only on chrM | `afquery query --db ./db/ --locus chrM:3243` | Mitochondrial coverage is easy to lose on its own: a BED naming it `chrMT` still matches on the autosomes, so nothing else looks wrong. Names are normalized (`MT`, `M`, `chrMT` → `chrM`) since v0.4.0 — databases built before that may need the capture BED re-checked. | | Technology filter | Remove `--tech` filter | Check if any samples match the requested technology. |