Skip to content

Unify the biallelic definition across GenotypeMatrix loaders and make biallelic restriction evident with warning - #164

Open
nspope wants to merge 2 commits into
nsp-multiallelic-trunkfrom
nsp-multiallelic-9
Open

Unify the biallelic definition across GenotypeMatrix loaders and make biallelic restriction evident with warning#164
nspope wants to merge 2 commits into
nsp-multiallelic-trunkfrom
nsp-multiallelic-9

Conversation

@nspope

@nspope nspope commented Jul 28, 2026

Copy link
Copy Markdown
Collaborator

Background

Two related problems. (1) Several entry points restrict to biallelic sites without telling the user,
so a biallelic-only statistic run on multiallelic data silently returns a partial result. (2) The
loaders disagreed on what "biallelic" even means, and on how to build a dosage -- so the same data
produced different genotype matrices depending on which loader you used, and none of that was
caught by the test suite because nothing compared the loaders.

This PR adds a shared BiallelicOnlyWarning, unifies the biallelic definition and the dosage rule
across all three GenotypeMatrix loaders (from_vcf, from_zarr eager, from_zarr streaming) and
from_haplotype_matrix, and adds a cross-loader parity harness.

The unified definition

  • Biallelic = at most two distinct present alleles (label / reference independent). A {0,2} site,
    or a reference-absent {1,2} site (possible after subsetting), is biallelic. This replaces the
    narrower "codes in {0,1}" notion that wrongly dropped genuinely-biallelic non-{0,1} sites.
  • Dosage = count of the chosen alt allele: allele 0 is reference, alt = highest present allele > 0,
    with a sentinel (out of range) when no alt is present so a monomorphic site counts to 0. This
    replaces summing allele indices (wrong for {0,2}) and gt > 0 (folds distinct alts; breaks on
    reference-absent {1,2}).
  • {0,1} input stays byte-identical, including a monomorphic-derived {1}-only site (dosage 2). This
    is what forced "allele 0 = reference" rather than "highest present allele" in the alt choice.
  • Both loaders retain monomorphic and all-missing sites (they load faithfully); dropping
    uninformative sites is the job of apply_biallelic_filter, not the loaders. The diploid stats
    tolerate such columns (grm excludes them via its polymorphic mask; pca/ibs guard every division).

The source of truth for biallelic sites is _biallelic_and_alt(ac) in genotype_matrix.py.

What changed

  • _warnings.py: BiallelicOnlyWarning(UserWarning) plus _warn_biallelic_only(n_dropped, *, context,
    stacklevel): coerces the count to int, no-ops at 0, emits a message leading with the context and
    dropped-site count.
  • GenotypeMatrix.from_haplotype_matrix: onto the shared definition. Counts the alt allele instead of
    summing indices; kepps {0,2}/{1,2}/monomorphic (previously {0,2}/{1,2} were dropped and index-sum
    gave impossible dosages like 2+2 -> 4); drops >= 3-allele sites + warns. {0,1} bit-identical.
  • GenotypeMatrix.from_vcf: moved off scikit-allel's is_biallelic_01() onto the shared definition.
    keeps {0,2}/{1,2}/monomorphic (is_biallelic_01 dropped all three); counts the alt allele; builds
    the allele-count matrix with numpy so scikit-allel's library footprint is now purely VCF parsing
    (read_vcf). Drops >= 3-allele sites + warns.
  • GenotypeMatrix.from_zarr (eager AND streaming, via build_genotype_matrix): onto the shared
    definition. Both zarr paths used a raw ploidy-sum with no classification, so a {0,2} site or any
    multiallelic site produced impossible dosages, silently, and disagreed with from_vcf. Now counts
    the alt allele. A >= 3-allele site has no valid 0/1/2 dosage, so its whole row is recoded to -1
    (present but fully missing) -- this preserves the row/chunk alignment the streaming path relies on
    (from_vcf drops the site instead). The warning is surfaced by the callers: the eager path with the
    exact count, the streaming path once per load on the first affected chunk.
  • admixture.patterson_d: warn with the count of >2-allele sites it already drops (num=den=0), once
    per top-level D computation (patterson_d / moving_patterson_d / average_patterson_d). D unchanged.
  • GenotypeMatrix.apply_biallelic_filter: docstring corrected. It does no allele classification (the
    matrix is biallelic by construction) -- it drops monomorphic (non-segregating) sites. Name kept
    for backwards compatibility; a rename is deferred as a follow-up.

Behaviour

  • {0,1} input is byte-identical across every loader, silent.
  • {0,2} / reference-absent {1,2} / monomorphic sites are now KEPT (and correctly dosed) by all
    loaders; previously from_vcf dropped them and the other paths mis-handled them.
  • A >= 3-allele site: from_vcf and from_haplotype_matrix drop it (fewer variants); from_zarr recodes
    its whole row to -1 (present but fully missing). All paths emit one BiallelicOnlyWarning.
  • All three GenotypeMatrix loaders now produce byte-identical dosage matrices from equivalent
    biallelic data, and stay consistent on multiallelic data.
  • The warning is silenceable by category (from pg_gpu import BiallelicOnlyWarning;
    warnings.filterwarnings("ignore", category=BiallelicOnlyWarning)).

Not in scope / deferred

  • LD / moments-LD paths keep the old "codes in {0,1}" notion by design (they are on their way to a
    separate multiallelic-correct LD rewrite): haplotype_matrix.apply_biallelic_filter (a moments
    is_biallelic_01 reimplementation), the moments-LD biallelic filter + alt construction, and the
    Rogers-Huff matrix dosage. These three call sites are documented in the LD subproject plan so the
    reconciliation is not lost. rogers_huff itself restricts to biallelic and will use this warning
    there.
  • apply_biallelic_filter is misnamed (it drops monomorphic sites, not a biallelic classification).
    Docstring fixed here; a rename or a clearer alias with soft-deprecation is a public-API change
    deferred as a follow-up issue.

Verification

  • tests/test_biallelic_only_warning.py: the class and helper (subclass, count/context message,
    numpy-scalar count, no-op at 0, silenceable); the shared _biallelic_and_alt definition and alt
    choice; from_vcf keeps {0,2}/{1,2}/monomorphic and drops the triallelic site with correct
    dosages; from_haplotype_matrix bit-identity on {0,1}, alt-code independence, reference-absent
    {1,2}, and the one-missing-haplotype -> fully-missing (-1) case; degenerate-column downstream
    safety (grm and pca_dosage stay finite on explicit monomorphic and all-missing columns);
    patterson_d warns once.
  • tests/test_loader_parity.py (new): drives from_vcf, from_zarr (eager), and from_zarr (streaming)
    from a single call_genotype block and asserts they agree byte-for-byte on biallelic data, stay
    consistent on multiallelic data (zarr recodes to a -1 row where from_vcf drops), and that the
    streaming warning fires exactly once across multiple chunks. This is the cross-loader coverage
    that was missing when the three loaders silently diverged.

@nspope
nspope changed the base branch from main to nsp-multiallelic-trunk July 28, 2026 03:34
@nspope
nspope requested a review from andrewkern July 28, 2026 03:34
nspope added 2 commits July 27, 2026 20:36
- Introduce a shared BiallelicOnlyWarning (UserWarning) and a _warn_biallelic_only
helper in _warnings.py, exported from the package. The warning marks a
statistical-domain restriction -- a biallelic-only statistic or loader dropping
multiallelic sites -- distinct from MultiallelicCapWarning (a kernel-capacity cap
on a multiallelic-capable statistic). The helper coerces the count to int, no-ops
at 0 so callers pass a raw count unconditionally, and leads the message with the
context and dropped-site count plus a silence recipe. Foundational for the sites
that restrict to biallelic (from_vcf, from_haplotype_matrix, patterson_d) and for
rogers_huff in the LD subproject; those emit sites follow in later commits.

- from_vcf filters to is_biallelic_01() sites and previously dropped everything else
silently. Emit BiallelicOnlyWarning with the dropped-site count at load, so a user
loading a multiallelic VCF sees how many sites were excluded. No change to which
sites are kept. Tests cover a VCF with a triallelic site (warns with count, loads
the biallelic subset) and an all-biallelic VCF (no warning).

- from_haplotype_matrix summed paired haplotype allele indices for every site, so a
site with an allele index >= 2 produced a bogus dosage (2+2 -> 4, or 1+2 -> 3, both
impossible). The conversion is documented as well-defined only for {0,1} alleles,
so this was silent garbage on multiallelic input. Drop sites whose max allele index
exceeds 1 (and their positions), and emit BiallelicOnlyWarning with the count.
Unlike the from_vcf and patterson_d warn sites -- which surface
filtering that already happened -- this ADDS a filter, a deliberate
behavior fix: multiallelic input now yields the biallelic subset plus a
warning instead of a matrix with impossible dosages.

- patterson_d (ABBA-BABA) is defined on biallelic SNPs and drops sites with >2
alleles across the four populations (num=den=0), previously silently. Emit
BiallelicOnlyWarning with the dropped-site count in _patterson_d_gpu, which is
called once per top-level computation (patterson_d / moving_patterson_d /
average_patterson_d all window the full num/den arrays afterward), so the warning
fires once per call, not per window. Warn-only: the D values are unchanged.
- Introduce genotype_matrix._biallelic_and_alt(ac): from an (n_var, K) allele-count
matrix it returns the biallelic mask (<= 2 distinct present alleles) and the
alt-allele index the 0/1/2 dosage should count (the highest-index present allele,
with an out-of-range K sentinel for monomorphic sites so the count is 0 and never
collides with the -1 missing code). This is the single source of truth for what
"biallelic" means and which allele the dosage counts -- label- and
reference-independent, so {0,2} and reference-absent {1,2} sites are handled
correctly. It is the per-individual form of patterson_d's alt-frequency choice.

- Move from_haplotype_matrix onto _biallelic_and_alt: a site is biallelic iff it has
<= 2 distinct present alleles (not codes in {0,1}), and the 0/1/2 dosage is the
per-individual count of the alt allele ((c1==alt)+(c2==alt), missing -> -1) rather
than a sum of allele indices. So a genuinely biallelic {0,2} site -- or a {1,2}
site with the reference dropped by subsetting -- is now KEPT with the correct
dosage instead of being dropped or turned into a bogus dosage; sites with >= 3
distinct alleles are dropped with the warning. {0,1} data is bit-identical.

- Refine _biallelic_and_alt so allele 0 is the reference and alt is the highest
present allele > 0 (sentinel only when no alt allele is present); this keeps {0,1}
bit-identical including a site fixed for allele 1, which the earlier
distinct-count sentinel would have zeroed.

- Move GenotypeMatrix.from_vcf onto the shared _biallelic_and_alt definition,
matching from_haplotype_matrix. That is: (A) Biallelic = at most two
distinct present alleles (was scikit-allel's is_biallelic_01). {0,2} and
reference-absent {1,2} sites are now KEPT and handled per-allele;
monomorphic and all-missing sites are retained too. Only sites with >= 3
alleles are dropped (with a BiallelicOnlyWarning).
(B) Dosage counts the chosen alt allele (highest present > 0) instead of summing
allele indices, so non-{0,1} codings are correct. (C) Build the
allele-count matrix with numpy; allel's library footprint is now purely
VCF parsing (read_vcf). The count ignores missing (-1) and matches
allel.count_alleles().

- build_genotype_matrix (shared by the eager and streaming zarr paths) computed a
raw maximum(gt,0).sum(axis=2) ploidy sum, so a {0,2} biallelic site or any
multiallelic site produced impossible dosages -- and it disagreed with the now
unified from_vcf loader. Now it classifies each site with the shared _biallelic_and_alt definition (at
most two distinct present alleles) and counts the chosen alt, matching from_vcf
and from_haplotype_matrix. A >= 3-allele site has no valid 0/1/2 dosage, so its
whole row is recoded to -1 (present but fully missing), which keeps the
row/chunk alignment the streaming path relies on (from_vcf drops the site
instead). A BiallelicOnlyWarning is surfaced by the callers: the eager path with
the exact count, the streaming path once per load on the first affected chunk.
{0,1} input (the common case) is byte-identical to before -- sum and counted-alt
agree there.

- tests/test_loader_parity.py is new: it drives from_vcf, from_zarr (eager), and
from_zarr (streaming) from a single call_genotype block and asserts they agree
on biallelic data and stay consistent on multiallelic data. This is the
cross-loader coverage that was missing when the three silently diverged.
@nspope
nspope force-pushed the nsp-multiallelic-9 branch from 4c63cd6 to 733e07b Compare July 28, 2026 03:40
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant