Multiallelic-correct windowed daf_hist / mu_sfs#160
Open
nspope wants to merge 7 commits into
Open
Conversation
The windowed daf_hist and mu_sfs derived their counts from dac = sum(max(hap, 0)), which sums allele indices (index-inflated on multiallelic sites) and used a fixed n_hap with a total-variant denominator (wrong under missing data). Rewrite the include fast path in windowed_statistics_fused to the per-allele convention of the tskit-pinned scalar diversity.daf_histogram / diversity.mu_sfs, so the windowed value over a window equals the scalar on that window's variants (windowed_statistics_fused_chunked delegates these stats to the same function). daf_hist now bins one frequency per present derived allele (count / n_valid) and normalizes per window; mu_sfs counts segregating and edge (site, allele) entries with per-site n_valid, fixing the fixed-n_hap misclassification that flagged a monomorphic-derived site with one missing haplotype as a near-fixed edge. A shared diversity._daf_bin_index makes the windowed scatter and the scalar histogram assign bins identically (diversity._histogram_from_dafs now bincounts through it).
The fused engine assigned each variant to a window with searchsorted on window ENDS, which places a variant sitting exactly on a boundary (pos == a window end == the next window's start) in the lower window -- disagreeing with win_start / win_stop, n_variants and the dedicated scatter engines, which are right-open [start, end). Assign instead to the last window whose start <= pos, with a pos < end containment check. Fixes the boundary misassignment for all the bin_idx-scattered fused stats (daf_hist, mu_sfs, mean_nsl) in one place.
Under missing_data='exclude', daf_hist and mu_sfs fell through to the WindowedAnalyzer fallback, whose registry did not know them, so they raised "Unknown statistic". Register them as a FREQ_STATS category that defers to the scalar diversity.daf_histogram / diversity.mu_sfs per window (which handle exclude and multiallelic correctly), so the windowed value equals the scalar over the same variants. daf_hist is vector-valued and expands into daf_bin_0 .. daf_bin_19 columns; an empty window yields 0.0 / an all-zero histogram, matching the scalar on an empty slice and the fused engine. No dispatcher change: exclude already routed to the fallback.
The windowed daf_hist bin count was hardcoded as 20 in two independent places (the fused engine's n_daf_bins and the fallback's n_bins), which had to agree or the include and exclude paths would emit different daf_bin_* columns. Define one module constant _DAF_N_BINS in windowed_analysis and use it in both engines; the tests reference the symbol rather than the literal. No behavior change.
daf_hist (normalized histogram) and mu_sfs (edge/segregating ratio) are not plain sums, so streaming correctness is not obvious. It holds because the streaming grid aligns windows to chunk boundaries (window_size divides align_bp; a window never straddles two chunks), so each window is computed whole within one chunk by the per-chunk eager path and the frames are concatenated -- no cross-chunk normalization. Add an equivalence test across chunk/window combinations (one-window-per-chunk through whole-fixture-single-chunk). No library change.
…e xfail) The streaming equivalence test used a no-missing fixture. Add a ~10%-missing vcz fixture and split the daf_hist/mu_sfs streaming test by missing_data mode: include matches eager (per-site n_valid is chunk-invariant), exclude is xfailed (strict=False). Streaming + exclude diverges from eager at chunk boundaries for ALL windowed stats -- plain pi shows the same start/n_variants disagreement -- a pre-existing streaming window-grid issue that was never exercised because the streaming tests only ran include.
Register mu_sfs in the cross-implementation parity matrix (scalar reference diversity.mu_sfs; fused and pyloop paths), so it is checked against the scalar across clean / missing-include / missing-exclude / multiallelic (the fused engine is auto-skipped under exclude, the fallback covers it). daf_hist is vector-valued and sits outside the scalar matrix, so add a dedicated test comparing its whole-region daf_bin_* vector to diversity.daf_histogram for the fallback (all modes) and the fused engine (include).
Collaborator
Author
|
The hardcoded number of bins for daf_hist for the streaming path isn't ideal ... but this is a niche statistic and I'm not sure of how to forward arguments correctly from windowed_analysis, so I'm leaving that alone |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
The windowed diploSHIC features daf_hist and mu_sfs derived their derived-allele counts from
dac = sum(max(hap, 0)), which sums allele indices, so on a multiallelic site the derived-allele
frequency was index-inflated and could exceed the sample size. They also diverged from the
tskit-pinned scalar diversity.daf_histogram / diversity.mu_sfs on plain biallelic data: a fixed
n_hap instead of per-site n_valid, raw counts instead of a normalized histogram, a total-variant
instead of a segregating denominator, and a missing-data edge misclassification in mu_sfs. And
exclude mode was unsupported for both (it raised "Unknown statistic").
This aligns both windowed features to the scalar functions so the windowed value over a window
equals the scalar computed on that window's variants, across biallelic / multiallelic and clean /
missing data. GenotypeMatrix input is biallelic by construction and out of scope.
What changed
Include (the fused fast path). windowed_statistics_fused (and its chunked twin, which delegates
these stats to it) now computes daf_hist and mu_sfs per allele via _memutil.allele_counts with
per-site n_valid, matching diversity: daf_hist bins one frequency per present derived allele
(count / n_valid) into (window, bin) and normalizes per window; mu_sfs counts segregating
(0 < count < n_valid) and edge (count == 1 or n_valid - 1) (site, allele) entries and returns
edge / segregating (0.0 when a window has none). A shared diversity._daf_bin_index makes the
windowed scatter and the scalar histogram assign bins identically.
Exclude. Previously unsupported (it fell through to the WindowedAnalyzer fallback, whose registry
did not know these stats). Registered as a FREQ_STATS category that defers to diversity.mu_sfs /
diversity.daf_histogram per window (which already handle exclude and multiallelic), with daf_hist
expanded into daf_bin_0 .. daf_bin_{_DAF_N_BINS - 1} columns. This keeps the package's standard
include-fused / exclude-fallback split, both pinned to the scalar.
Window boundary. While here, the fused engine's per-variant window assignment (bin_idx) was
changed from searchsorted on window ends (which placed a variant sitting exactly on a boundary in
the lower window) to the right-open [start, end) assignment used by win_start / win_stop,
n_variants and the scatter engines. This fixes the boundary misassignment for all bin_idx-scattered
fused stats at once (daf_hist, mu_sfs, mean_nsl).
Constant. The windowed daf_hist bin count, hardcoded as 20 in two places, is now the single module
constant _DAF_N_BINS. Ideally this should be settable but that's fiddly.
Behaviour change
These are now the per-allele definitions of the two features: multiallelic-correct, using
allele_counts and per-site n_valid, the same per-allele convention as the SFS, max_daf and the
scalar diversity.daf_histogram / mu_sfs. Because the features now follow that one consistent
per-allele convention rather than their old index-summed / fixed-n_hap / raw-count one, their
output changes on all HaplotypeMatrix input, biallelic included -- a deliberate redefinition of the
two diploSHIC feature vectors. daf_hist becomes a normalized per-window histogram over
present derived alleles (no monomorphic zero-frequency spike), and mu_sfs uses the segregating
denominator with correct per-site edge classification, returning 0.0 rather than NaN for a window
with no segregating entries. "exclude" mode is newly supported (was an error).
Verification
diversity.daf_histogram / diversity.mu_sfs (pinned to tskit's polarised allele frequency spectrum)
are the oracle; the contract is windowed == scalar over a window. Covered: whole-region windowed ==
scalar for biallelic / multiallelic and clean / include-missing (fused) and exclude-missing
(fallback), including a monomorphic-derived-plus-missing site that the old fixed-n_hap edge test
misclassified; a multi-window test with variants on interior boundaries (the boundary fix); a
streaming == eager equivalence across chunk/window combinations (include, and include with missing
data); and both stats added to the cross-implementation parity harness (mu_sfs in the scalar
matrix, daf_hist as a dedicated vector check) across clean / include / exclude / multiallelic.
Deferred
a pre-existing streaming window-grid issue never exercised before (streaming tests only ran
include). This is logged as Streaming + missing_data='exclude' diverges from eager at chunk boundaries (all windowed stats) #159