Anchor every oracle's effect null on the regions its assay measures - #163
Merged
Conversation
…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>
…an 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>
…ually 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>
#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>
…n 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>
…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>
…e 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>
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>
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>
…nded 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>
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.
A variant-effect percentile answers "how unusual is this effect, compared to what?". For a localised assay, uniformly random genomic positions are the wrong answer — almost none carry the signal the model reads, so the pseudocount damps their log-ratios toward zero, the null's body collapses, and real effects pile up at the top of the scale where the column no longer discriminates.
The concrete symptom: at SORT1, half of Enformer's chromatin-accessibility rows had a percentile pinned at exactly 1.0000. Now zero.
What the fleet looked like
Measured by reading which pass consumes each position source — not assumed:
snps = random + dhsTwo of those last three already sampled cCREs for their baseline pass. An asymmetry inside a single oracle is harder to defend than any difference between oracles. Now 8 of 8.
Result — per-track tail ratio, new null vs old
Calibration held. Strong TSS-proximal liver eQTLs (GTEx v8, tissue-matched): AlphaGenome RNA 0.781 → 0.778, CAGE 0.659 → 0.625. Both in band, both 0% saturated.
It is a union at doubled N, not a mixture — that difference is the finding
The first attempt held the position count fixed and gave cCRE 25% of it. It made things worse. The statistic that decides whether a percentile still discriminates is the null maximum, and a maximum grows with the number of draws — so splitting a fixed budget shortens every component's tail:
Enformer TF saturation went 25% → 92%. Keeping each component at full size makes the union's maximum exactly
max(max_gene, max_cCRE), so it is provably never worse than the better component for any layer — a guarantee, not a measurement. The gene-anchored half reproduces the previously shipped counts exactly (1,200 / 1,200 / 1,980 / 720), so the cCRE half is purely additive.LegNet got a different set on purpose. It's a 200 bp promoter MPRA model with
window_bp=None, so the sampled position is the whole thing modelled. The cCRE catalogue is 62% dELS against 2% PLS, and DHS summits track accessibility not promoter identity — either would give a promoter model a null made mostly of enhancers. So TSS ±250 bp 40%, PLS 30%, pELS 15%, uniform 15%.A silent bug this surfaced
Enformer wrote its internal
spec_key(DNASE,CHIP_TF, …) as its per-row layer while AlphaGenome wrote canonical names. A composition keyed onchromatin_accessibilitymatched 472 of AlphaGenome's rows and 0 of Enformer's 5,313 — silently, for the one oracle where the fix had been measured to help. Same class as #122 and #144. Fixed with a sharedcanonical_layer()that raises rather than falling back, all builders routed through it, plus a guard test.Also #144 instance five:
discovery.pyhand-rolled centred windows at two sites, summing 4–5 bins against a null built over 3, in the path that ranks cell types. The umbrella was closed on a hand-enumerated list of four; the new guard searches instead of trusting an enumeration.Three measurements I got wrong first
median(new)againstmedian(old)— the median track of each set, which needn't be the same track. Per-track it is 1.38×.All three are the same error: a summary statistic that didn't mean what I wanted. Acceptance is now per-track p99/p99.9. The maximum is still the right thing to report — it's literally "pinned at 1.0" — but the wrong thing to tune on.
Not fixed, and not claimed to be
histone_marks(20%) and Enformertf_binding(25%) keep whatever their better component gives. Both need a per-track population — that mark's own domains, that factor's own peaks — a different design, not a different fraction.Verification
*_countsconsecutive-integer run (Enformer tracks are not ranked against the same variant set (effect_counts vary 9600–9606) #123's fingerprint).Nothing published to HuggingFace.
🤖 Generated with Claude Code