Skip to content

regen(examples): make the blog's two worked examples reproducible + correct the draft - #136

Merged
lucapinello merged 3 commits into
mainfrom
regen/2026-08-01-with-cebp-full-track-set
Aug 1, 2026
Merged

regen(examples): make the blog's two worked examples reproducible + correct the draft#136
lucapinello merged 3 commits into
mainfrom
regen/2026-08-01-with-cebp-full-track-set

Conversation

@lucapinello

Copy link
Copy Markdown
Contributor

Closes the Track B work: every quantitative claim in the blog post is now either reproducible from a committed artefact, or corrected, or explicitly marked as unverifiable.

Analysis B (rs9504151 / CDYL) rested on nothing committed

Every number traced to one markdown table in audits/2026-06-16_blogpost_reproduction_report.md, describing a GPU run whose inputs and outputs were never saved. Only 9 rsIDs from that locus appear anywhere in git history — not the 56-variant proxy set — and the AlphaGenome track list behind "rank #1, composite 0.995" was recorded nowhere. The fine-map path had not run since 2026-06-17.

Both inputs are now committed (ld_proxies.tsv, assay_ids.txt), so --only cdyl reproduces it in ~2 min with no LDlink token. And it does reproduce:

audit recorded re-run today
composite 0.995 0.991
alt effect −1.363 −1.362
Rank Variant Composite Largest effect
1 rs9504151 1.00 0.991 −1.362 · DNASE
2 rs658325 1.00 0.686 −1.109 · DNASE
3 rs386522231 0.87 0.652 −1.035 · DNASE
4 rs62384944 0.93 0.441 +0.293 · DNASE

So the audit's numbers were right — they just couldn't be checked by anyone.

Resolves the 54-vs-56 discrepancy the audit left open (it reported "54 (CEU, SNVs)" and "rank #1/56" in adjacent cells, never reconciled): snvs_only=True56 (1 sentinel + 55 proxies), snvs_only=False64. Same set, counted with and without the sentinel.

Four things the README now says that the post does not: the ranking is a composite of raw values not percentiles; rs62384944 sits in a cluster of 11 variants at that same r², so its rank is more arbitrary than implied; chorus computes no PIP (0.51 is from the Sniff preprint); and LegNet cannot adjudicate this locus (no lung fibroblast in K562/HepG2/WTC11), so the post's "LegNet ranked second" has no supporting run anywhere.

Analysis A: two of the blog's four C/EBP figures were uncheckable

validation/SORT1_rs12740374_with_CEBP was a byte-for-byte duplicate of variant_analysis/SORT1_rs12740374, and its README claimed RNA+/- tracks it never scored. Widened to 11 tracks:

blog committed now
CEBPB +3.1 +3.044
CEBPA +2.8 +2.764
CEBPG +2.3 +2.269
CEBPD +1.9 +1.818

CEBPG and CEBPD previously appeared in no example at all. RNA at +0.718 = 1.65-fold makes the "well under 2-fold" claim checkable too — and the README notes that AlphaGenome's RNA-seq background sits almost entirely below |log2FC| = 0.1, so read the log2FC there, not the percentile.

The blog's drift explanation is measurably wrong

The post attributes a +1.24 → +1.37 gap to "floating-point arithmetic differences across GPU/CPU runs":

hypothesis measured
run-to-run float noise 0.0000% — three identical runs, +1.375621 each, spread 0.000000e+00
expexpm1 (#113) −0.12%, opposite direction
1 bp → full window +0.318 → +1.376 (4.33×) — verified at 4ad7be7

ChromBPNet is bit-exactly deterministic. (AlphaGenome/JAX is not — but this claim is about ChromBPNet.) +1.24 is the original article's figure, predates the audit, and its provenance is unknown; the draft now says that rather than substituting a tidier story.

Draft corrections (7 more)

audits/2026-006_chorus_blogpost.revised.md items 9–15: oracle count seven → eight models / nine registered with the missing Cherimoya row added; tools 23 → 24; AlphaGenome resolution qualified (2,915 of 5,731 are 1 bp, all 2,816 CHIP are 128 bp); Sei 21,907 vs 40 scoreable; background description corrected (not "thousands" uniformly, and not gnomAD); display axis clips at 3× p99.

A note at the top records that the published post diverged from this file — it shipped the pre-review Enformer window while taking the post-review tool count.

The size guard earned its keep

The CDYL report renders at 25.70 MB, above the 20 MiB ceiling — and tests/test_committed_examples.py caught it before it shipped, the same failure that cost a rejected push earlier this week.

I did not raise the ceiling. The cause is a real gap: #129's budget is per track; nothing caps a report total, so 21 tracks × 2 alleles is ~42× a single-track panel. Filed as #135 with a proposed report-level budget; until then the generator skips the HTML unless CHORUS_WRITE_LARGE_HTML=1, so it cannot silently recreate a rejected artefact.

626 passed, 8 skipped, 0 failed · zero ref-allele mismatch warnings across all regenerated examples

🤖 Generated with Claude Code

lucapinello and others added 3 commits August 1, 2026 17:36
…ect the blog draft

TWO THINGS, both about making a public claim checkable against a committed
artefact.

1. validation/SORT1_rs12740374_with_CEBP was a byte-for-byte duplicate of
   variant_analysis/SORT1_rs12740374 — same position, alleles, gene and
   track list — and its README claimed tracks it did not score:

     "Tracks used: DNASE, CEBPA ChIP, CEBPB ChIP, H3K27ac, CAGE+/-, RNA+/-"

   Neither RNA strand was ever scored, and only two of the four C/EBP family
   members were. Widened to 11 tracks (adds CEBPG, CEBPD and both HepG2 polyA
   RNA strands), which makes the README true rather than editing the claim
   away, and gives the example a purpose distinct from its twin. 63 -> 123
   scored rows.

   The whole family responds, which is the point of a validation case:

     CHIP:CEBPB:HepG2   +3.044      blog says +3.1
     CHIP:CEBPA:HepG2   +2.764      blog says +2.8
     CHIP:CEBPG:HepG2   +2.269      blog says +2.3
     CHIP:CEBPD:HepG2   +1.818      blog says +1.9
     DNASE:HepG2        +1.332
     CHIP:H3K27ac:HepG2 +1.258
     ATAC:HepG2         +0.732
     RNA:HepG2          +0.718  (best of 58 rows)

   All four C/EBP figures the blog cites now reproduce within rounding, from
   a committed file — previously CEBPG and CEBPD appeared in no example at
   all, so two of them were uncheckable.

   RNA +0.718 is a 1.65-fold predicted increase: real, and well short of the
   >12-fold measured SORT1 mRNA difference between genotypes. The README now
   says so explicitly, and warns that AlphaGenome's RNA-seq effect background
   sits almost entirely below |log2FC| = 0.1, so the RNA *percentile* carries
   almost no information even where the effect does. Read the log2FC column.

2. audits/2026-006_chorus_blogpost.revised.md — seven further corrections,
   each re-measured against the code (items 9-15 in its changelog):

     - "seven oracles" -> eight models / nine registered names, and Cherimoya
       was missing from the spec table entirely despite landing 2026-07-30
       with 1,518 calibrated tracks.
     - MCP tools 23 -> 24 (score_ism landed after the first revision).
     - AlphaGenome "1 bp resolution" qualified: 2,915 of 5,731 tracks are
       1 bp; all 2,816 CHIP tracks are 128 bp. And 5,168 of 5,731 have
       percentile backgrounds — AUDIT_CHECKLIST.md:205 requires always saying
       which number is meant.
     - Sei "21,907 chromatin profiles" qualified: chorus aggregates to 40
       scoreable classes, 0.18% of that head.
     - The background description was wrong twice over: it is not "thousands
       of random variants" uniformly (1,697-1,909 for AlphaGenome vs 18,672
       for ChromBPNet) and it is not gnomAD — no code samples gnomAD.
     - Display axis: 1.0 is the per-bin p99 but the axis runs to 3.0, so the
       old phrasing implied 1.0 was the ceiling.

   And the Analysis-A reproducibility note is rewritten, because the
   published post's explanation is measurably wrong. It attributes a
   +1.24 -> +1.37 gap to "floating-point arithmetic differences across
   GPU/CPU runs". Measured:

     - ChromBPNet is BIT-EXACT: three identical runs of DNASE:HepG2 at
       chr1:109274968 all gave +1.375621, spread 0.000000e+00. There is no
       float noise to attribute anything to. (AlphaGenome/JAX *is*
       nondeterministic, but this claim is about ChromBPNet/TensorFlow.)
     - exp -> expm1 accounts for -0.12% and in the opposite direction:
       +1.373940 vs +1.375621 from one model run with only that transform
       swapped. It exactly explains the audit's +1.374 becoming today's
       +1.3756, and nothing more.
     - The 1 bp auto-region bug is a different comparison: +0.318 vs +1.376
       (4.33x). Verified by running at 4ad7be7, which gives +0.317674 and
       reproduces the audit's recorded +0.32.

   So +1.24 is the ORIGINAL ARTICLE's figure, it predates the audit, and its
   provenance is unknown — most likely a different window. The draft now says
   that instead of blaming float noise, which also wrongly implied the number
   could not be trusted.

   A note at the top of the draft records that the *published* post diverged
   from this file — it shipped the pre-review Enformer window (196,608 bp,
   actually 393,216) while picking up the post-review tool count. Reconcile
   against the file, not the published text.

Zero ref-allele mismatch warnings across all four regenerated examples.

622 passed, 7 skipped, 0 failed.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
… CDYL)

Every quantitative claim in the blog post's Analysis B traced to one markdown
table in audits/2026-06-16_blogpost_reproduction_report.md, describing a GPU
run whose inputs and outputs were never saved. Across the whole git history
only 9 rsIDs mention rs9504151's locus — not the 56-variant proxy set the
numbers were computed from — and the AlphaGenome track list behind "rank #1,
composite 0.995, alt effect -1.363" was never recorded at all. The fine-map
path has not been exercised since 2026-06-17.

Committed, so the claim becomes checkable:

  examples/walkthroughs/causal_prioritization/CDYL_rs9504151/
    ld_proxies.tsv   56 rows, the LDlink LDproxy response with the exact
                     query in its header
    assay_ids.txt    the 21 AlphaGenome lung-fibroblast tracks (CL:0002553)

This is the same trick SORT1_locus already uses — it inlines its 11 proxies as
LDVariant literals, which is precisely why SORT1 regenerates with no LDlink
token and rs9504151 could not. 56 literals would be unreadable inline, so this
is a data file plus a `_load_ld_fixture` reader, and the file records the query
parameters rather than leaving them to be inferred.

RESOLVES THE 54-vs-56 DISCREPANCY the audit left open. It reported "54 (CEU,
SNVs)" in one cell and "rank #1/56" in an adjacent one and never reconciled
them. Measured against LDlink today:

    snvs_only=True   -> 56   (1 sentinel + 55 proxies)
    snvs_only=False  -> 64

so the two figures are the same set counted with and without the sentinel,
modulo a small 1000G panel drift since June. `snvs_only` is now recorded
explicitly in the fixture header rather than inferred. r2 spans 0.8020-1.0000
with 27 variants at exactly 1.0.

Also confirms the blog's "rs62384944 (r2~0.93)" — it is present. Worth noting
for the write-up that it sits in a cluster of **11** variants at essentially
that same r2, so singling out its rank is more arbitrary than the prose
implies.

scripts/regenerate_remaining_examples.py gains `regen_cdyl_finemap` and
`--only cdyl`, wired into the `all` and `skip_discovery` task lists so this
example cannot become an orphan the way three others did.

NOT YET RUN. AlphaGenome is a gated HuggingFace model and the locally stored
token was just rotated, so the model load fails with "Invalid user token".
The fine-map re-run needs `hf auth login --force` first. Expect the recorded
ChromBPNet numbers (composite 0.896, effect -0.985) to move when it does —
four corrections landed after that audit (expm1 counts, CDF regen, percentile
denominator, BPNet CHIP) and both `max_effect` and `ref_activity` read off the
corrected counts.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
rs9504151 ranks #1 of 56, composite 0.991, largest effect -1.362 on
DNASE:fibroblast of lung. Against the figures
audits/2026-06-16_blogpost_reproduction_report.md recorded but never made
checkable:

    composite    0.995  ->  0.991
    alt effect  -1.363  -> -1.362

So the audit's AlphaGenome numbers were right. They simply could not be
verified by anyone, because the 56-variant proxy set and the 21-track
assay_ids list were never saved. Both are committed now, so
`--only cdyl` reproduces this with no LDlink token in ~2 minutes.

    Rank  Variant        r2    Composite  Largest effect
    1     rs9504151 *    1.00  0.991      -1.362 DNASE
    2     rs658325       1.00  0.686      -1.109 DNASE
    3     rs386522231    0.87  0.652      -1.035 DNASE
    4     rs62384944     0.93  0.441      +0.293 DNASE
    5     rs9504169      0.83  0.421      +0.005 RNA

rs62384944 at rank 4 matches the audit's AlphaGenome figure exactly (the
blog's "rank 4" was graded PARTIAL only because ChromBPNet put it at 5).

The README records four things a reader needs and the post does not say:

  - the ranking is a composite of RAW values, not percentiles
    (causal.py:756-788, weights 0.35/0.25/0.20/0.20), so #83's percentile
    problems do not enter it — but neither does any magnitude calibration;
  - rs62384944 is not uniquely at r2 ~ 0.93. It sits in a cluster of 11
    variants at essentially that r2, so singling out its rank is more
    arbitrary than the prose implies;
  - chorus computes no posterior inclusion probability at all. PIP = 0.51
    is from the Sniff preprint (`grep -rni "posterior_inclusion|susie|
    polyfun" chorus/` -> zero hits);
  - LegNet cannot adjudicate this locus — its panel is K562/HepG2/WTC11,
    no lung fibroblast — so the published post's "LegNet ranked second"
    has no supporting run anywhere, and the repo's own corrected draft
    explicitly declines to use LegNet here.

The ChromBPNet IMR-90 figures (composite 0.896, effect -0.985) are NOT
re-run and are flagged stale: four corrections landed after that audit and
both max_effect and ref_activity read off the corrected counts.

NO COMMITTED HTML, and the guard is why. The causal report renders at
25.70 MB, above the 20 MiB ceiling tests/test_committed_examples.py
enforces — the guard added yesterday caught it before it shipped, which is
the same failure mode that cost a rejected push earlier this week. I did
not raise the ceiling: its headroom over today's largest artefact
(14.7 MiB) is what makes it useful.

The cause is a real gap rather than a bad example. #129's budget is PER
TRACK (4,000); nothing caps a report total, so 21 lung-fibroblast tracks x
2 alleles is ~42x a single-track panel:

    per-track budget    4,000
    tracks x alleles    21 x 2 = 42
    worst case          168,000 features

Filed as #135 with a proposed report-level budget. Until then the
generator skips the HTML unless CHORUS_WRITE_LARGE_HTML=1, so it cannot
silently recreate an artefact the guard rejects. The JSON/MD/TSV carry
every number; only the browser panel is missing.

626 passed, 8 skipped, 0 failed.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@lucapinello
lucapinello merged commit b6193ee into main Aug 1, 2026
1 check passed
@lucapinello
lucapinello deleted the regen/2026-08-01-with-cebp-full-track-set branch August 1, 2026 18:35
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