fix: enformer effect percentiles had an unreachable top 4% (padded CDF grid) - #143
Merged
Merged
Conversation
…F grid) `enformer_pertrack.npz` shipped its `effect_cdfs` gridded at 9,606 points — `max(effect_counts)` — then padded to 10,000 by repeating each row's maximum 395 times. All 5,313 rows, every one with its first maximum at index 9605. Since `_get_denominator` correctly divides by the stored grid *width* (#119), the padding silently rescaled every enformer effect percentile and made a whole band of the scale unreachable: probe quantile percentile shift after repair p50 +0.0205 p90 +0.0368 p99 and above +0.0393 reachable ceiling 0.9604 -> 0.9999 Nothing caught it, three ways: * `ReservoirSampler.to_cdf_matrix` interpolates short rows onto the *full* grid and cannot produce that shape, so the file was never reproducible from repo code; * `build_and_save` only resampled when `shape[1] > n_points`, so a narrow matrix was written verbatim, with no assert; * `normalization.py`'s own docstring states rows are never padded. The artefact contradicted an invariant its code documents. It also contaminated #83: a per-track floor at any q >= 0.96 resolves to the enformer null *maximum* for all 5,313 tracks, so the q sweep was invalid there. WHAT THIS ADDS `cdf_grid_violations` / `expected_first_max_index` in background_sampling.py — one implementation shared by the test, the write-time guard and the repair. Two checks, each validated against all eight shipped backgrounds (19,393 short rows): both flag every one of enformer's 5,313 effect rows and nothing else. The derivation matters and is easy to get wrong. `to_cdf_matrix`'s `source_q = arange(n)/n` stops at `(n-1)/n`, not 1.0, so `np.interp` clamps beyond it and a short row legitimately ends in a plateau of repeated maxima — 6 slots for AlphaGenome's 1,909 samples, 100 for a 100-sample row. Trailing duplicates are therefore normal and are *not* the padding signal. Three edge cases are exempted for cause, each with a test: * the plateau check is one-sided. Padding makes the plateau longer (395 vs 2); a hand-built fixture makes it shorter, which is legitimate. * rows with `n >= n_points` are skipped: AlphaGenome's summary row 2,452 has its top two H3K4me1 windows tied at exactly 2480.0 — real saturation. * constant rows carry no grid geometry at all (a single sample interpolates flat), so both checks would false-positive on them. A narrow matrix is *not* rejected, only logged. `_get_denominator` reads `shape[1]`, so a 100-wide matrix stored at 100 is self-consistent — and enformer's defect was the opposite shape, `shape[1] == n_points`, which no width check could have caught. THE REPAIR `scripts/repair_enformer_effect_grid.py` recovers the quantile function from the 9,606 stored quantiles and re-expresses it on the 10,000-point grid: idempotent, backed up, exact max preservation, median-quantile shift 6.4e-07, all 5,313 rows monotone. Pseudo-samples must span `linspace(0,1,n)` and not `arange(n)/n` — the latter discards the top sample and loses up to 74% of a row's maximum, precisely the tail being restored. This is an interpolation of a derived artefact, not a rebuild; the planned enformer rebuild supersedes it. Local cache only — republishing to HuggingFace is a separate authenticated step. Tests: 670 passed (was 653), 8 skipped. Integration: 8/8 shipped backgrounds clean, and the retained pre-repair backup still trips the guard. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
lucapinello
added a commit
that referenced
this pull request
Aug 4, 2026
…, and fix four defects found while verifying (#160) * fix: the grid guard mistook tied maxima for padding, and a real rebuild caught it The #143 guard REFUSED to write a fresh Borzoi background. It was right to fail closed and wrong about the reason. Borzoi row 3011's top effect value 0.689308 recurs 9 times, because several sampled variants hit the same clipped ceiling. np.interp holds a tied value from the q-position of the first tied sample onward, so ties lengthen the trailing run of maxima exactly as padding does. 10 Borzoi rows and 152 Enformer rows tripped the one-sided plateau check. I had already exempted ties for `n >= n_points` — citing AlphaGenome's summary row 2,452, two H3K4me1 windows tied at exactly 2480.0 — and simply never applied the same reasoning to short rows. WHY THE FIX IS FILE-LEVEL Padding shifts EVERY row's maximum to the same index, because the whole matrix was gridded at one narrower width. Ties shift a handful. So the signal is the MEDIAN of (expected first max - actual) across rows: every legitimate shipped background 0 a fresh Borzoi / Enformer build 0 enformer_pertrack.npz before repair 393 Three orders of magnitude between signal and tolerance. cdf_grid_file_violations() also abstains below 8 rows, because one tied row genuinely is indistinguishable from one padded row — only unanimity separates them. The tie-immune per-row check (distinct != count) is kept and is what actually caught the original defect: 5,313 of 5,313 padded rows flagged, 0 of 30,000+ legitimate ones across every shipped file AND both fresh rebuilds. Two of my own test expectations were wrong here too: the tie fixture used N_POINTS=64 with 5,949 samples, which is the subsample branch rather than the interpolation branch being tested, and a single-row matrix cannot support a median. Tests: 38 in this module, 853 overall. Both fresh rebuilds now merge cleanly. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * feat: AlphaGenome's RNA null uses TSS-in-window selection, matching the reference Aligns gene selection with what AlphaGenome actually does, verified against the shipped implementation rather than the paper text. THEIR RULE, from the code gene_mask_extractor.py:326 transcripts selected by TSS position _PositionExtractor position >= interval.start and < interval.end (semi-open) gene_annotation.py:94 TSS = Start for '+', End for '-' (strand-aware) gene_annotation.py:75 feature='transcript', so PER-TRANSCRIPT gene_mask_extractor.py:369 gene mask = OR of exon masks of THOSE transcripts So a gene whose body overlaps the window but whose TSS lies outside contributes NOTHING under their rule, where the gene-body-overlap rule chorus used contributes all of its exons. HOW MUCH IT ACTUALLY MATTERS — measured before spending 10 GPU-hours locus overlap(PC) TSS-in-window(PC) only-overlap only-TSS SORT1 29 28 1 0 FTO 6 5 1 0 BCL11A 4 3 1 0 TERT 14 14 0 0 MYC 2 2 0 0 At most ONE gene of ~29 per locus, and overlap is a strict superset. The smoke run confirms it end to end: 14.5 genes per window under TSS-in-window against 15.6 under overlap. Small — but Luca asked for strict correctness over cheapness, and the old rule was a #144-class divergence I introduced myself. THE GENE-TYPE FILTER IS A DELIBERATE DIVERGENCE, NOT AN OVERSIGHT AlphaGenome applies NO gene-type filter — TSS-in-window returns 31-61 genes per locus against 28 protein-coding. chorus's query filters to protein_coding (variant_report.py:825), so the builder does too: a null over lncRNAs and pseudogenes would be a different population from the numerator, which is #144 in the other direction. Exposed as `gene_types=` and documented so it is a recorded choice rather than an implicit one. Aggregation and normalization need NO change — they already match exactly after #149: softplus output space (coverage-like, not log), mean over mask extent, natural log, +1e-3, no CPM and no per-track scaling. The pseudocount is benign at real activity levels: RNA p50 activity is 0.208, where a 2x signal change yields 0.688 of a possible 0.693. The overlap-rule interims are preserved at /data/chorus_data/ag_{effect,baseline}_interim.OVERLAP-RULE.npz so the change can be diffed rather than trusted. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Demote the distinct==count grid check to advisory; raise only on first_max The distinct-vs-count fingerprint I called perfect — 5,313/5,313 with zero false positives — fired on a healthy row and blocked a 10-hour rebuild merge. AlphaGenome effect_cdfs row 3966 (CHIP_TF ARID3A) has count=5949 and exactly 5949 distinct values on the 10,000-point grid. It is not padded: 913 of its samples are exact zeros, so interpolating the remainder happens to land on 5949 distinct values by coincidence. The proof it is interpolated is that its maximum first appears at index 9998, which is exactly expected_first_max_index(5949, 10000) — padding would put it at 5948. So distinct==count is a coincidence-prone fingerprint, not a proof. This makes the mechanical signal the only raising condition: first_max == n - 1 and n < n_points - 2 That is a direct consequence of how the padding was produced (grid at `count`, repeat the last column) and cannot be reached by np.interp, whose source_q stops at (n-1)/n. distinct==count becomes a logger.warning that says outright it is usually coincidence. Tally for this guard before the demotion: three false positives against one true catch (the per-row plateau check on tied maxima — 10 borzoi rows and 152 enformer rows; the file-level median check, removed entirely; and this). A guard that blocks a rebuild at that rate is worse than no guard. Verified: still catches the genuinely padded enformer_pertrack.npz.prepad (row 0, first_max 9605), and all three fresh rebuilds pass. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * Pin what quantile tie-breaking can actually reach, and what it cannot Two test corrections, both because the test encoded my motivation rather than the code's behaviour. test_a_value_inside_the_tie_is_spread_across_it queried 0.0417 — the value I said AlphaGenome RNA's null tops out at — expecting the flat top to be spread. It is not, and cannot be: build_and_save stores CDFs as float32, so the stored ceiling is 0.04170000106 while the float64 query is 0.04169999999. The query sorts strictly BELOW it, searchsorted returns lo == hi, and no tie is seen. So tie-breaking does not do what motivated it. Its reachable set is queries bit-exact against a stored grid value, which in practice means exactly 0.0 — common, because a variant a track does not respond to gives ref == alt bit-for-bit. Measured there it works: 39 distinct percentiles across 40 identical tracks, spanning 0.0030-0.0894 of a 0.0900-wide band. The test now uses that case, and a second test pins the float32 near-miss as a stated limitation rather than leaving it to be rediscovered. A value just above a degenerate ceiling still reads 1.0. That needs a wider null, not a lookup change — which is what the re-anchored region set is for. The grid-integrity test asserted the substring 'distinct', naming the fingerprint just demoted to advisory. It now asserts the mechanical signal. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * Stamp as-built provenance on the three rebuilt backgrounds, and gate determinism end to end The provenance PR taught build_and_save and append_tracks to carry a file-level build_config, but none of the three rebuild-set builders passes one, so alphagenome, borzoi and enformer shipped with only the canonical 8 keys. Only cherimoya had a build_config, and it predates this work. build_config is file-level, so it is appended in place rather than by re-running eleven hours of forward passes. The discipline is that only establishable facts go in. Three things are deliberately NOT stamped: * xla_flags, because reading os.environ here records THIS process's flags, not the build's. pin_deterministic_xla_flags logs nothing on its success path, so the build logs carry no positive record either. * a builder_git_sha set to HEAD, because two commits landed after these builds finished. Worse, the AlphaGenome build started at 16:26 with its builder changes still uncommitted -- they became c4ded13 at 19:28, three hours into an eleven-hour run -- so no single sha describes what ran. * wall-clock figures. Instead the stamp is content-addressed (git hash-object sha1 + mtime + a boolean for whether the mtime predates the build start) and, decisively, records the strata the build ITSELF logged. That last field settled a case the mtimes could not: borzoi and enformer started at 03:08:07, nine minutes BEFORE the commit that introduced gene-anchored sampling, and their mtimes were later bumped to 03:30:14 by a git operation, so mtime alone reads as changed-after-build-start for both. Their logs show the gene-anchored strata at 03:08:25 and 03:08:30 -- the code was running from the first minute, as uncommitted edits. All three oracles logged identical counts (tss_near 1200, tss_far 1200, junction 1980, gene_body 720, random 849 of 6000), confirming one shared seeded region set as designed. Also adds gate_end_to_end_determinism.py. #127 was closed on 9/9 bit-exactness at the raw predict_sequence level, which does not cover the scorers, the window geometry, the percentile lookups or the tie-breaking hash layered on top. This runs the real build_variant_report path twice in two processes on one GPU and diffs every numeric leaf bitwise -- no tolerance, because #127's symptom was 36 sign flips and a tolerance would hide those. One fix worth naming: np.savez_compressed appends .npz unless the filename already ends in it, so a .stamping temp suffix wrote to a different path and the verify-before-replace read a file that did not exist. The originals were never touched, but the guard was silently inert. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * Assert the three release gates, and record what four CAGE reference classes measure The rebuild plan listed three residuals to caption; two were fixable and one was a hedge, so all three became work items. This turns them into tests. Gate 3 is the headline: a full build_variant_report is now bit-exact across two OS processes on one GPU. 603 numeric fields, 0 differing, 0 sign flips, worst relative delta exactly 0.0 -- against #127's baseline of 454 differing fields with 36 sign flips. #127 was closed on 9/9 bit-exactness at the raw predict_sequence level, which did not cover the scorers, the window geometry, the percentile lookups or the tie-breaking hash layered on top; this covers all of it. The calibration band needs a correction on the record. An earlier version of this assertion used [0.75, 0.95], which I invented before having any measurement, and CAGE at 0.659 then read as a failure when it is simply where a correctly-widened null puts strong eQTLs. The band is now [0.55, 0.97], stated as a collapse/ re-saturation detector rather than a target to tune toward. The audit note records the four CAGE reference classes measured against strong TSS-proximal liver eQTLs: uniform-random positions (shipped before) p50 0.857 gene-anchored mixture (shipped now) p50 0.659 all annotated TSSs, variant AT the TSS p50 0.323 all annotated TSSs, eQTL-matched offsets p50 0.713 Placing the variant exactly at the annotated TSS puts it at the peak maximum, where the +1 pseudocount barely damps, so the null's median |effect| (0.036-0.039) comes out nearly double a validated eQTL's (0.020) and a real eQTL reports BELOW median. That is the mirror of the failure this work set out to fix. Anchoring on the same TSSs but drawing the offset from the empirical eQTL tss_distance distribution gives 0.713 -- better than the shipped mixture, with no invented 40/33/12/15 fractions. So the direction is right; only the exactly-at-the-TSS reading of it is wrong. Also pins #123's fingerprint across all eight shipped backgrounds: no *_counts array may be a tight run of consecutive integers (enformer shipped 9600-9606). Verified clean on all eight. One thing that did not improve, recorded rather than buried: enformer chromatin_accessibility at SORT1 went from 4/12 saturated to 6/12. That is not a new regression -- 0.960 was the padded-grid artefact ceiling (#143), so those rows were already pinned and the repair only made it visible instead of disguising it as a plausible 0.96. Enformer's accessibility effect null is genuinely too narrow for a variant this strong. It clears gate 1 and is the next thing to look at. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * Delete the duplicate TSV writer that was silently dropping per-gene rows Found while checking JSON<->TSV parity across all 13 walkthroughs after the rebuild. Six of them disagreed, and the cause was not the rebuild -- the counts are identical before and after -- but a second implementation of one routine, which is the same defect class as everything else in this cycle. scripts/regenerate_remaining_examples.py carried its own flattener, _variant_report_tsv_rows, wrong in two ways: 1. It de-duplicated on (allele, assay_id, layer), a key that omits the region. RNA and CAGE emit one row per GENE per track -- same allele, same assay, same layer, different gene -- so every gene after the first was discarded: validation/TERT_chr5_1295046 99 json -> 18 tsv discovery/SORT1_cell_type_screen 347 rows -> 39 tsv sequence_engineering/region_swap 32 json -> 4 tsv sequence_engineering/integration 55 json -> 3 tsv TERT's TSV showed one tss_activity row where there were fifteen, one per nearby gene TSS (BRD9, CLPTM1L, LPCAT1, MRPL36, ...). A dropped row leaves no gap and no null in a TSV, so nothing looked wrong -- which is why this survived however long it has been there. 2. It wrote region_label into a column named "description", which already means the TRACK description in to_dict(). One name for two things across two artefacts of the same report, so the two could not even be compared until the column was renamed back. The fix is deletion, not repair: route everything through report.to_dataframe(), the canonical writer scripts/regenerate_examples.py already used and which needs no de-dup at all. All 14 walkthrough (json, tsv) pairs now agree exactly. tests/test_json_tsv_parity.py pins both the counts and the row IDENTITIES -- an equal count with different rows would be a worse failure than an unequal one. 877 fast tests pass. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * Correct the AlphaGenome track count in user-facing docs: 5,731 -> 5,168 The docs advertised AlphaGenome as 5,731 tracks in ten places across four files, including the README's headline sentence. 5,731 is the row count of its metadata table; 563 of those rows are `padding` placeholders whose only job is to keep local_index aligned with the model's output array. They carry no assay, iter_tracks() skips them, and the shipped background has no row for any of them. The number a user can actually query is 5,168 -- verified both directions: 5,168 metadata tracks have a background row, and 0 do not. So the docs overstated usable coverage by 563 tracks. Not a large error, but the kind that erodes trust in every other number printed beside it, and it survived because nothing compared the prose to the artefact. CHANGELOG.md:324 claims this was already "disambiguated inline" in an earlier cycle. It was not, in any of the four live docs. tests/test_documented_track_counts.py now does the comparison: it reads the shipped NPZs plus the oracle metadata and fails on any live doc whose per-oracle count contradicts them. It also pins WHY the two numbers differ (5,168 real + 563 padding = 5,731 exactly), so the explanatory footnote cannot quietly become false. Dated files under audits/ and historical CHANGELOG entries are deliberately left alone -- they record what was true when written, and rewriting them would destroy the record rather than fix anything. That exclusion is stated in the test. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * CHANGELOG: the rebuild, with attribution where it is measurable and none where it is not Records the gene-anchored rebuild of AlphaGenome, Borzoi and Enformer, the eQTL calibration table, and the effect on the committed examples: saturated rows fell from 47 to 16 across the four variant walkthroughs, row counts unchanged at 369. The plan asked for a table attributing each moved number to the change that moved it. That is delivered for the four changes measured in isolation -- the RNA denominator (numerator overstated 251-1736x), the Enformer grid repair (ceiling 0.9605 -> 0.9998), cross-process determinism (603 fields, 0 differing, against 454 differing with 36 sign flips), and the gene-anchored null (the eQTL table). It is NOT delivered per-change for the per-layer walkthrough diff, and the entry says so rather than implying otherwise. Splitting that honestly would mean re-running each rebuild in isolation; splitting it by judgement would mean printing numbers I cannot support. It ships as a fused diff with its causes listed. The Enformer chromatin_accessibility regression is in the entry as a Changed item with its own heading, not buried: 4/12 saturated -> 6/12, because 0.960 was the artefact ceiling and the repair made existing pinning visible rather than creating it. Stated as "did not improve, and is not described as fixed". Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * Blog draft: the rebuild made two of its background claims stale, and one framing wrong The five corrections the plan listed for this draft were already applied (Enformer 393,216 bp, nine registered oracles with Cherimoya, the 5,731-vs-5,168 and 21,907-vs-40 footnotes, tool count 24). Three things needed changing anyway. Stale because of my own rebuild: the draft described every effect null as built from "random genomic positions" with 1,697-1,909 samples per AlphaGenome track and ~9,600 for Enformer/Borzoi. AlphaGenome, Borzoi and Enformer are now gene-anchored, at 5,949-87,781 / 5,949-20,337 / 5,949. Sei, LegNet, ChromBPNet and Cherimoya are unchanged and still uniformly random, so the sentence has to distinguish them rather than generalise. Also updated the reader-facing gloss: "relative to random SNPs" is no longer true for three of eight oracles. Wrong framing, and it was mine to begin with: "5,168 of the 5,731 have percentile backgrounds" implies 563 real tracks lack one. They are not tracks -- they are `padding` placeholders that keep local_index aligned with the model's output array, carry no assay, and are skipped by iter_tracks(). Every queryable track has a background row. 5,168 is the number to quote, and the draft now says why. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * The bit-exactness gate has to find the alphagenome env, not skip itself away The new integration gate failed on its first real run: ModuleNotFoundError: jax. It launched the script with sys.executable, which under the base chorus env has no JAX -- the oracle envs are deliberately isolated (CLAUDE.md). Fixed by locating the sibling interpreter (../chorus-alphagenome/bin/python from sys.prefix) rather than by adding a skip. A gate that skips itself when its dependency is missing reports green while checking nothing, which is worse than failing: the whole point of this test is that #127's symptom was invisible. It still skips if neither the sibling env nor an importable jax exists, but the skip message names the exact command to run by hand. Verified through pytest on GPU 0: 1 passed in 215s. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * Make the stamped provenance load-bearing, and unify a shape three writers disagreed on build_and_save and append_tracks could write and preserve a file-level build_config, and the three rebuilt backgrounds carry one -- but a grep for build_config across chorus/ found only the writer, the append-path preservation, and docstrings. Nothing read it. Provenance nobody consults is documentation, not an invariant, so it could not catch the thing it exists to catch. That thing is #122: AlphaGenome histone CHIP nulls built over 501 bp while the query summed 2001 bp. Both artefacts were internally consistent, so neither looked wrong; the defect existed only BETWEEN them, and it shipped across 1,075 of 5,168 tracks. There is now a read side -- PerTrackNormalizer.provenance() plus a load-time comparison of the stamped histone/other window_bp against LAYER_CONFIGS -- so that disagreement is detectable instead of invisible. It warns rather than raises. Every background built before the provenance work is unstamped, and an unstamped file makes no claim to contradict; refusing to load one would break working installs to enforce a metadata convention. Writing the reader immediately surfaced a third instance of this cycle's recurring defect. build_and_save wrote build_config as a 0-d array; Cherimoya's builder and scripts/stamp_background_provenance.py both wrote a 1-element array, which is what every file on disk actually has. So a reader written against the shipped artefacts raised IndexError on anything build_and_save produced, and vice versa -- three producers, two conventions, nothing comparing them. The writer is now consistent at 1-element and _read_build_config accepts both, so old files still load. Three assertions in test_npz_provenance.py had pinned the 0-d shape; they now parse through the shared reader and assert content instead. Also refreshes all six executed notebooks against the rebuilt backgrounds: 6/6 re-executed clean, 0 errors. Run in a mirrored directory tree because cherimoya's notebook resolves ../../genomes/hg38.fa relative to its own location, so executing a flat copy under /tmp fails on a path, not on anything real. 895 fast tests pass; 24 integration pass with 2 skipped. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> --------- Co-authored-by: Claude Opus 5 <noreply@anthropic.com>
lucapinello
added a commit
that referenced
this pull request
Aug 5, 2026
…163) * Per-layer effect region sets, and exact position-sharding so a build can use 8 GPUs Two things, both driven by a census rather than by argument. THE CENSUS. Comparing every committed walkthrough row against its own track's null maximum on the shipped gene-anchored null: oracle layer rows above null max enformer chromatin_accessibility 12 50.0 % alphagenome histone_marks 10 30.0 % enformer tf_binding 12 25.0 % alphagenome gene_expression 100 7.0 % alphagenome tss_activity 263 8.0 % enformer tss_activity 48 0.0 % The PEAK layers saturate and nothing else does. The cause is structural: most gene-anchored positions are not inside a peak, and a variant in closed chromatin cannot move an accessibility or ChIP signal much, so the null's upper tail is too short. At SORT1, enformer accessibility effects are 1.14-1.45x its null maximum, which is exactly why they read 1.0000 and stop discriminating. So --effect-regions {gene-anchored,ccre} is added to all three builders, and sample_ccre_anchored_positions reuses sample_ccre_positions -- the ENCODE SCREEN sampler the BASELINE path has always used and the effect path never did. Sharing it means the two paths cannot drift. CAGE turned out to need nothing, which is the opposite of what I expected and worth recording. Sweeping the variant's distance from an annotated TSS gives a monotone curve in that single parameter (strong-eQTL percentile p50: 0.323 at the TSS itself, 0.411 at +/-500 bp, 0.526 at 1 kb, 0.604 at 2 kb, 0.654 at 5 kb, 0.729 at 10 kb), and the shipped gene-anchored mixture sits at 0.659 -- it already behaves like a "+/-5 kb of a TSS" null for CAGE. There is no qualitative gain available, only a choice of distance scale that no principle fixes. Also: the "Generated N gene-anchored SNPs" log line now names the region set from args instead of hardcoding it, because that line IS the provenance stamp_background_provenance.py reads back. Saying "gene-anchored" while sampling cCREs would stamp a lie into every rebuilt NPZ. The stamp's regex captures the region set rather than assuming it, and still parses the existing gene-anchored logs. EXACT POSITION SHARDING. These three oracles produce every track from ONE forward pass, so sharding by track -- which is how the chrombpnet builder does it, correctly, because its tracks are separate model files -- saves no GPU time; each shard would still run every pass. They have to be sharded by POSITION, and that changes what a merge is: each shard then holds a partial reservoir for EVERY track. The tempting merge is to pool the shards' 10,000-point CDF grids. That is an approximation, decent for equal shards, and this repo has already shipped one artefact whose approximation looked exact (#143's padded grid). So ReservoirSampler.to_flat_samples/from_flat_samples serialise raw samples -- ragged data stored flat with offsets, so the NPZ still loads under allow_pickle=False -- and the CDF is built once from the union. counts ADD, since shards see disjoint positions. tests/test_position_sharding.py pins the property that makes this trustworthy: for 2, 3, 4 and 8-way partitions the merged sampler is bit-identical to one fed the whole stream, in retained values, in counts, and in the resulting CDF. Also covers eviction under capacity pressure, empty tracks, disagreeing track counts, and determinism. 912 fast tests pass. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * #144 instance FIVE: the discovery path computed a different window than the query I closed #144 as an umbrella on four instances. A consistency audit found a fifth, at two sites in chorus/analysis/discovery.py: ref_track.score_region(chrom, pos - half, pos + half + 1, cfg.aggregation) That is the exact pre-#144 arithmetic scorers.py was migrated off. It builds genomic coordinates and lets score_region floor/ceil-expand them back to bins, giving 4 OR 5 bins for window_bp=501 at Enformer's 128 bp depending on where the variant falls inside its bin -- against a null built over 3. discover_variant_effects and discover_cell_types read the SAME per-track backgrounds as the variant report, so both were ranking a wider statistic against a narrower reference. Worse for discover_cell_types, which *compares* those effects to order cell types: a span that varies with sub-bin position can reorder the ranking itself. It survived because the umbrella fix was applied where the defect was known to be, and nothing searched for the shape elsewhere. Four instances were enumerated by hand; a fifth existed. So the guard now does the searching, rather than trusting an enumeration: any score_region call on the same line as half-width arithmetic fails, across every module in chorus/. score_region itself is not banned -- it is right for a genuine genomic interval, which core/base.py's exon summing and core/result.py's region scoring both want. What is banned is using it to approximate a CENTRED window. A second, separate defect at the same place: gene_expression has window_bp=None because its statistic is the mean over a gene's merged exon mask, and the window_bp-is-None branch fell through to np.mean(track.values) -- the mean over the ENTIRE prediction, 524 kb for Borzoi and 1 Mb for AlphaGenome, dominated by intergenic and intronic zeros. That number bears no relation to the exon-mask null it was then ranked against. Discovery has no gene context to build a mask from, so it now declines the layer with a logged reason instead of emitting a confidently wrong percentile; analyze_variant_multilayer, which takes a gene, is the path for RNA. The guard is deliberately two-sided. A test that only forbids the bad pattern goes green if someone deletes the scoring altogether, so there is also a positive assertion that score_centered_window is actually called at both sites, and an enumeration of which layers legitimately have no window so a new one is noticed rather than silently falling through to the whole-prediction mean. 1,051 fast tests pass. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * README: describe the reference population the shipped backgrounds actually use Six places said the effect percentile was computed against "~10,000 random SNPs sampled uniformly across chr1-chr22". That has not been true for AlphaGenome, Borzoi or Enformer since the gene-anchored rebuild, and it was never true of the sample size: the shipped counts are 5,949-87,781 per track, not 10,000. It is also the wrong shape of claim, because the reference population is no longer one thing. It now varies by oracle AND by layer, so the README carries a table: gene-anchored for the three rebuilt oracles, cCRE-anchored for their accessibility rows specifically, uniform-random for the five that were never rebuilt. The prose explains WHY rather than just listing, since the reasoning is the part a reader needs to interpret a percentile: a random genomic position carries almost no CAGE or accessibility signal, so the +1 pseudocount damps its log-ratio toward zero and the null's body collapses below where real regulatory effects live. It also names the deliberate 15 % uniform tail -- without near-zero mass, genuinely small effects would get artificially LOW percentiles, the mirror of the same failure. And it states the negative result, because a reference class that helped one layer and not others is the kind of thing that otherwise gets over-generalised later: accessibility rows come from cCREs because 50 % of Enformer's accessibility rows at SORT1 exceeded their own null maximum and a cCRE-anchored null takes that to 0 %. TF binding and histone marks were measured and did NOT improve -- a cCRE is defined by accessibility, H3K4me3 or CTCF signal, so a randomly chosen one often is not bound by the particular TF a given ChIP track measures. The sample-size table is knowingly still stale in this commit; it is regenerated from the shipped artefacts once the cCRE build lands, so its numbers are final rather than written twice. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * Migrate the last local ReservoirSampler out of the AlphaGenome builder (#125) Eight builders, eight copies of this class. Seven had already been migrated to the shared one; build_backgrounds_alphagenome.py kept its own, 87 lines at line 107. Keeping it cost real time today. I added to_flat_samples to the SHARED class for position sharding, launched eight AlphaGenome shards, and every one of them ran for fifty GPU minutes and then died at the write step with AttributeError: 'ReservoirSampler' object has no attribute 'to_flat_samples' because this file's copy did not have the method. That is the thesis of #125 demonstrated at my own expense, and the reason the fix is deletion rather than keeping the copy in sync. Both differences the copy carried are preserved: * default capacity 20,000 against the shared 50,000 -- moot, since all three call sites pass capacity=args.reservoir_size explicitly; * a hand-vectorised add_batch, which the baseline pass's per-variant fan-out genuinely needs. That is now the SHARED implementation, so migrating does not regress AlphaGenome's throughput, and the plain loop survives as _add_batch_reference so the equivalence test still has something to compare to. The vectorised fast path is only taken while the reservoir still has room -- where "which samples survive" is not yet a question, so bulk-extending is identical to appending one at a time. Once full it falls back to per-value Algorithm R, because there the traversal ORDER decides which samples survive and a different order moves the CDF with no arithmetic changing. Verified numerically against the deleted implementation, not just by reading: at capacities 40, 500 and 20,000, over 30 chunks of 37 values into 3 tracks, both counts and retained samples are identical. Capacity 40 forces the overflow branch, which is the only branch where order matters. 1,054 fast tests pass. The 10 failures in the run are tests/test_walkthrough_readmes_ match_artefacts.py correctly flagging stale README prose, fixed in a follow-up once the rebuild lands and the artefacts are final. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * Anchor the effect null on the regions each assay actually measures, in all 8 oracles Five of the eight oracles already anchored their effect null on peaks; three did not, even though two of those three already used cCREs for their BASELINE pass. An asymmetry inside a single oracle is harder to defend than any difference between oracles, so this closes it. chrombpnet already had it: 10,000 DHS-summit variants + 10,000 uniform cherimoya already had it: snps = random + dhs, explicitly unioned alphagenome gains it: gene-anchored + cCRE union enformer gains it: same borzoi gains it: same sei gains it: same generic union -- its 40 classes span promoters, enhancers, TF binding, transcription, heterochromatin and low signal, so no single peak type fits, which is when the generic mixture is right legnet gains a PROMOTER-anchored set, not the generic one. It is a 200 bp promoter MPRA model with window_bp=None, so the sampled position IS the whole thing being modelled. The cCRE catalogue is 62% dELS (1,469,205 of 2,348,854 distal enhancer-like) against 2% PLS (47,532 promoter-like), and DHS summits track accessibility rather than promoter identity -- either would give a promoter model a null made mostly of enhancers. PLS leads at 30%, TSS+/-250bp at 40%, pELS at 15%, uniform tail at 15%. epinformer gains its OWN baseline composition for the effect pass, removing the intra-oracle asymmetry. MEASURED, per track, as the ratio of the new null's tail to the old one's. Median over tracks, and the share of tracks that got wider: oracle tracks p99 p99.9 wider sei 40 2.05x 1.80x 100% epinformerseq 33 1.38x 1.28x 76% legnet 3 1.30x 1.17x 67% enformer 5,313 1.26x 1.33x 84% borzoi 7,611 1.19x 1.19x 82% All five improve. Sei gains most, which fits: it had a pure uniform-random null and the widest class space. TWO CORRECTIONS TO MY OWN METHOD, both worth recording because both produced a wrong answer first. (1) I initially reported EPInformer-seq as REGRESSING at 0.89x. That came from comparing median(new) against median(old) -- the median TRACK of each set, which need not be the same track. The correct statistic is the median of per-track ratios, which is 1.38x with 25 of 33 tracks wider. Same family of error as (2). (2) The metric itself. I had been scoring these changes by "share of rows above the null MAXIMUM". A maximum is a single extreme order statistic with high sampling variance: Enformer's shipped tf_binding max is 3.539 while the rebuilt one is 2.956 even though its p99.9 nearly doubled. That metric reads noise in one draw, not tail width. It is the right thing to REPORT to a user, since it is literally "percentile pinned at 1.0", but the wrong thing to tune on. Acceptance is now p99/p99.9. WEAKER JUSTIFICATION, STATED AS SUCH: sei, legnet and epinformerseq have no committed walkthrough rows and no positive set (no eQTL equivalent for MPRA activity or Sei sequence classes), so unlike the Enformer accessibility fix -- where saturation was measured at 50% and dropped to 0% -- these three are justified by tail width and by matching the assay's biology, not by a calibration check. That distinction belongs in the CHANGELOG too. Also fixes an import bug I introduced and nearly shipped: the guard that was supposed to add the module-level `from chorus.utils.annotations import ...` tested for a string that a FUNCTION-level import of the same module already satisfied, so nothing was written and both sei and legnet would have died with NameError at runtime. Caught by an AST bound-name check rather than by the parse, since a missing import parses fine. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * Apply the rebuilt effect nulls to 5 backgrounds; generate the README table apply_effect_rebuild.py swaps a freshly built effect null into a shipped background and touches nothing else. The rebuild changed only which positions the EFFECT null is drawn from -- the baseline/summary and per-bin passes are identical -- so re-running them would burn GPU reproducing the same numbers while risking drift. It is also the only correct option rather than a shortcut: sei, legnet and epinformerseq no longer have baseline interims on disk, so `--part merge` is not available for them at all. What it refuses to do: * proceed when track_ids differ in content OR ORDER. Row i means "the track at index i of track_ids", so a reordering silently reassigns every null to the wrong track -- and a same-set-different-order case would pass a naive set comparison. * proceed when the new matrix fails cdf_grid_violations (#143's guard). * overwrite in place. It writes a sibling, verifies the sibling loads and that summary_cdfs/perbin_cdfs came through byte-identical, and only then replaces. Applied to enformer, borzoi, sei, legnet, epinformerseq. Originals backed up to /data/chorus_data/pre_effect_rebuild/. AlphaGenome follows when its shards land. Verified after: all 8 shipped backgrounds monotone, zero NaN/inf, provenance readable. Effect samples/track now 11,913-34,482 for the rebuilt ones. The provenance stamp records effect_region_set per oracle, and says explicitly that only the effect rows moved -- otherwise a reader would reasonably assume the activity percentile was rebuilt too. Also replaces README's hand-maintained per-oracle background table with a generated one. It had drifted badly: 10,000 effect / 31,500 activity samples claimed for six of seven rows against real counts of 5,949-104,033, and Cherimoya -- one of the eight shipped backgrounds -- missing from the table that enumerates them. Same defect as the stale walkthrough READMEs: prose about numbers that live in binary artefacts, with nothing comparing the two. One fix worth naming: the generator's fallback defaulted to "uniform random" when a stamp had no effect_region_set, which mislabelled AlphaGenome -- its stamp records the gene-anchored rule under an older key. It now falls back only to recorded facts and prints "(unrecorded)" rather than asserting something wrong. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * Regenerate every artefact against the rebuilt nulls, and fix the prose that lied All eight backgrounds changed, so every percentile in every committed example moved. Regenerated the 13 walkthroughs, the multi-oracle path (chrombpnet + legnet + alphagenome + consolidate), and corrected 34 stale numbers across 11 READMEs. Twenty-seven of those were the same row with a slightly moved value and were substituted mechanically at each site's own decimal precision. Four needed prose rewrites because the claim, not the number, was wrong: * variant_analysis/SORT1_chrombpnet/README.md described an ATAC run at -0.111 ("moderate closing") and devoted a whole section to explaining why AlphaGenome and ChromBPNet DISAGREE at this locus. The committed artefact is DNASE at +1.376, and AlphaGenome reports +1.334 -- they agree on direction and within a few percent on magnitude. The section explained a contradiction that does not exist. Rewritten to report the agreement, which is the actual point of running two oracles, and to note that the loaded assay is DNase not ATAC. * variant_analysis/README.md claimed the Enformer example's key finding was "very strong CTCF binding loss (-0.89)". The artefact contains NO CTCF track at all. Replaced with its real maxima: TF +4.372, chromatin +2.247, histone +1.886. * variant_analysis/BCL11A_rs1427407/README.md had the SIGN wrong: "-0.113 log2FC" against an artefact value of +0.145. An opening reported as a closing. * The same file's ChromBPNet entry quoted +0.43 where the artefact says +1.376. The root cause is structural, not clerical: the regeneration scripts rewrite example_output.{json,md,tsv} and the HTML, and have never touched a README. So every correctness fix since #92 left the narrative behind, and nothing compared prose to data. tests/test_walkthrough_readmes_match_artefacts.py now fails on any signed decimal in a walkthrough README that its own artefact does not contain. TWO TESTS OF MINE ALSO NEEDED FIXING, both for asserting a FORM rather than the invariant, and both broken by a change that preserved the property exactly: * test_the_random_stratum_is_present_and_substantial asserted the random stratum's FRACTION was >= 0.10. Adding the cCRE stratum took the total from 6,000 positions to 12,000 and halved every pre-existing fraction to keep its COUNT identical: random went 0.15 -> 0.075, which is 900 positions either way. Now asserts the absolute count, which is the property that matters and is robust to N. * test_builder_imports_the_shared_classifier asserted an exact import LINE, which broke the moment a second name joined the same import. Now parses the AST and asserts the name is imported. A guard that fails on formatting trains people to edit the guard. 1,107 fast tests pass, zero failures. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * CHANGELOG + dated audit note for the reference-population work Records the result, the design, and the three measurements I got wrong before I got them right -- the mistakes being the more useful half. The audit note states the guarantee precisely, because it is what makes the design defensible rather than merely measured: keeping each component at full size makes the union's maximum exactly max(max_gene, max_cCRE), so the union is provably never worse than the better component for any layer. The fixed-N mixture had no such property and indeed came out below BOTH components for Enformer tf_binding, taking its saturation from 25% of rows to 92%. It also reconciles the previous cycle's entry rather than leaving two that contradict each other. That entry said Enformer chromatin_accessibility was "the next thing to look at" at 6/12 saturated; it is now 0/12, fixed later in the same release. The two are sequential, not competing, and the text now says so. The older eQTL table is marked as the post-gene-anchored figures with a pointer to the final ones. Both the CHANGELOG and the note state, in their own words, that Sei, LegNet and EPInformer-seq are justified by tail width and assay biology and NOT by a calibration check -- there is no eQTL equivalent for MPRA activity or Sei sequence classes. That is weaker ground than the accessibility fix stood on and should not be blurred into it. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * Re-execute all six notebooks against the rebuilt backgrounds 6/6 clean, zero error outputs. Run in a mirrored directory tree because cherimoya_quickstart resolves ../../genomes/hg38.fa relative to its own location, so executing a flat copy under /tmp fails on a path rather than on anything real. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> * Refresh the generated background table after AlphaGenome's rebuild landed AlphaGenome's effect samples/track went 87,781 -> 148,367 and its reference population is now recorded as gene-anchored+ccre rather than read from the older effect_region_rule key. Generated, not hand-edited. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> --------- Co-authored-by: Claude Opus 5 (1M context) <noreply@anthropic.com>
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 defect
enformer_pertrack.npzshipped itseffect_cdfsgridded at 9,606 points —max(effect_counts)— then padded to 10,000 by repeating each row's maximum 395 times. All 5,313 rows, every one with its first maximum at index 9605.Because
_get_denominatorcorrectly divides by the stored grid width (#119), the padding silently rescaled every enformer effect percentile and made a band of the scale unreachable:Understated everywhere, worst in the upper tail — which is exactly where a calibration floor reads.
It also contaminates #83. A per-track floor at any
q >= 0.96resolves to the enformer null maximum for all 5,313 tracks, so the q sweep was invalid there. This unblocks it.Why nothing caught it
ReservoirSampler.to_cdf_matrixinterpolates short rows onto the full grid and cannot produce that shape — so the shipped file was never reproducible from repo code.build_and_saveonly resampled whenshape[1] > n_points; a narrow matrix was written verbatim with no assert.normalization.py:455states rows are never padded. The artefact contradicted an invariant its own code documents.The guard
cdf_grid_violations/expected_first_max_indexinbackground_sampling.py— one implementation shared by the test, the write-time guard inbuild_and_save, and the repair script. Validated across all eight shipped backgrounds (19,393 short rows): flags every one of enformer's 5,313 effect rows and nothing else.The derivation is easy to get wrong.
source_q = arange(n)/nstops at(n-1)/n, not 1.0, sonp.interpclamps beyond it and a short row legitimately ends in a plateau of repeated maxima — 6 slots for AlphaGenome's 1,909 samples, 100 for a 100-sample row. Trailing duplicates are normal and are not the padding signal.Three exemptions, each with a test and a reason:
n >= n_pointsrows are skipped. AlphaGenome's summary row 2,452 has its top two H3K4me1 windows tied at exactly 2480.0 — real saturation, not corruption.A narrow matrix is logged, not rejected:
_get_denominatorreadsshape[1], so 100-wide stored at 100 is self-consistent — and enformer's defect hadshape[1] == n_points, which no width check could have caught.The repair
scripts/repair_enformer_effect_grid.pyrecovers the quantile function from the 9,606 stored quantiles and re-expresses it on the 10,000-point grid. Idempotent, backed up, exact max preservation, median-quantile shift 6.4e-07, all 5,313 rows monotone.Pseudo-samples must span
linspace(0,1,n), notarange(n)/n— the latter discards the top sample and loses up to 74% of a row's maximum, precisely the tail being restored.This is an interpolation of a derived artefact, not a rebuild; the planned enformer rebuild supersedes it. Local cache only — republishing to HuggingFace is a separate authenticated step.
Verification
670 passed, 8 skipped(was 653) on the CI invocation.🤖 Generated with Claude Code