Figure 3: Polygenic risk rescales loading; carrier variants reshape composition¶
Runs locally against the authors' results (CODE). With RUN = False the cells show the current local figures; set RUN = True to regenerate each panel from the script named. No participant data are stored in this notebook; keep outputs cleared before committing (python ../tools/check_no_ids.py).
Panel a Polygenic score rescales signature loading¶
Script: src/reparam_fig3/fig3a_prs_theta_share_by_backbone.py
Environment: BACKBONE=slope3way
Upstream: triplet training checkpoints
Script header (inputs, recipe, outputs)
GENERATED by _make_fig3a_share.py -- SHARE-OF-DISEASE-LOADING variant (see that docstring).
Fig 3a -- polygenic risk RESCALES loading: mean theta on a disease's signature over age, top vs bottom 1% of
that disease's PRS. Pairs as in the published panel: T1D->S15, CAD->S5, AF->S0, breast->S19 (women only).
Written 2026-09-28 because the generator of the published Fig 3a could not be located in either repo. It is
run on the PUBLISHED backbone first as a check (its values should reproduce the published panel), then on the
TRIPLET backbone.
BACKBONE=delta0 theta = softmax(lambda_nokappa_v3_delta0_400k.pt) -- the published full-fit lambda
BACKBONE=slope3way theta reconstructed batch-by-batch from results_slope3way (own per-batch gamma + delta):
lambda_k = r_k(t) + zb(G37) g_lvl_k + t*zb(G37) g_slope_k + delta_k ; healthy row = shifted
a=0.6 disease-free curve + zb(G37) g_health + delta_H ; t = arange(T)/(T-1); zb = batch z-score.
(the reconstruction used for fig_carrier_stacks_by_backbone.py / the coherence rebuild)
PRS strata use the RAW standardised score (G_matrix column), not the PC-residualised one the triplet trains on,
so the SAME people form the top/bottom 1% on both backbones.
INPUTS data_for_running/{G_matrix,baselinagefamh_withpcs,reference_trajectories}; prs_names.csv; backbone files
OUTPUT reparam_fig3/fig3a_prs_theta_<BACKBONE>.{png,pdf} + _values.csv
RUN BACKBONE=delta0|slope3way /opt/miniconda3/envs/new_env_pyro2/bin/python claudefile/reparam_fig3/fig3a_prs_theta_by_backbone.py
PHI group means over >=4,000 people.
Panel b Carrier effects reshape signature composition (12 genes, incl. S8 fits)¶
Script: src/reparam_fig3/make_monogenic_reshape_atlas_triplet.py
Environment: SWEEPDIR={CODE}/reparam_fig3/clinvar_sweep_gp_slope3way_fixed, REGROW_MLH1={CODE}/reparam_fig3/fig_colon_backbone_slope3way_sweepfmt_slope3way_fF_etat.npz, FLAGSHIPS=LDLR,APOB,TTN,MLH1,MSH6,MSH6_S8,BRCA1_S8,BRCA2,BRCA2_S8,SERPINA1,COL4A4,COL4A3, MAXPANELS=12, OUTSUFFIX=_slope3way_fixed_shortlabels, SHORTLABELS={CODE}/reparam_fig3/fig3b_short_labels.csv
Upstream: atlas: src/reparam_fig3/run_atlas_triplet_fixed.sh (fit_clinvar_sweep_gp.py); S8 fits: run_brca1_s8_triplet.sh and the BRCA2/MSH6 S8 fits
Script header (inputs, recipe, outputs)
MONOGENIC side of the dissociation: eta(t) heatmaps (carrier-vs-non-carrier expression over age).
Reads the time-varying atlas (clinvar_sweep_gp/*_etat.npz, produced by fit_clinvar_sweep_gp.py) and renders:
(A) FLAGSHIP per-gene heatmaps: the scoped signature's selected diseases (rows) x age (cols), colored by
eta(t) -- the carrier expression offset. eta>0 means the carrier state makes the OCCUPIED signature
express that disease MORE at that age. This is the "signature expression reshapes over time" visual:
the cascade lights up with gene-specific TIMING (LDLR: chol->CAD->MI; BRCA: breast+cascade; TTN: HF/AF).
eta(t) IS the carrier-vs-non-carrier contrast (phi_tilde_carrier - phi_tilde_control = eta on scoped cells).
(B) TIMING ATLAS: across ALL genes, the home-disease peak-eta age (when the gene acts) -- one strip.
Contrast make_polygenic_scale_contrast.py: PRS only SCALES loading (level shift), never adds this age-structure.
INPUT : clinvar_sweep_gp/*_etat.npz (eta_t (D,T), gate (D,), inc_excess (D,T), Dpen (D,), peak_age, ages, ...)
OUTPUT: fig_monogenic_reshape_flagships.{png,pdf} + fig_monogenic_timing_atlas.{png,pdf} (+prov).
RUN via new_env_pyro2. Aggregates only (fitted eta, no EIDs).
TRIPLET COPY (2026-10-07) of make_monogenic_reshape_atlas.py (that file and ~/mendelnoulli-paper/analysis/ are left
untouched). run_reshape_atlas_triplet.sh passes SWEEPDIR / OUTSUFFIX / FLAGSHIPS / REGROW_MLH1, which the original
does NOT read -- run with them it silently read the bar sweep and overwrote the bar figures (restored 2026-10-07).
This copy reads all four: SWEEPDIR = atlas dir; OUTSUFFIX appended to every output name (required when SWEEPDIR is
set, so bar figures are never overwritten); FLAGSHIPS = comma list, PINNED (no back-fill); REGROW_MLH1 = registry
eta npz for the MLH1 row.
Panel c Peak age of the carrier effect on each home disease¶
Script: src/reparam_fig3/make_monogenic_reshape_atlas_triplet.py
Environment: SWEEPDIR={CODE}/reparam_fig3/clinvar_sweep_gp_slope3way_fixed, REGROW_MLH1={CODE}/reparam_fig3/fig_colon_backbone_slope3way_sweepfmt_slope3way_fF_etat.npz, FLAGSHIPS=LDLR,APOB,TTN,MLH1,MSH6,MSH6_S8,BRCA1_S8,BRCA2,BRCA2_S8,SERPINA1,COL4A4,COL4A3, MAXPANELS=12, OUTSUFFIX=_slope3way_fixed_clean, TIMING_CLEAN=1
Upstream: same atlas as panel b; one fit per gene
Script header (inputs, recipe, outputs)
MONOGENIC side of the dissociation: eta(t) heatmaps (carrier-vs-non-carrier expression over age).
Reads the time-varying atlas (clinvar_sweep_gp/*_etat.npz, produced by fit_clinvar_sweep_gp.py) and renders:
(A) FLAGSHIP per-gene heatmaps: the scoped signature's selected diseases (rows) x age (cols), colored by
eta(t) -- the carrier expression offset. eta>0 means the carrier state makes the OCCUPIED signature
express that disease MORE at that age. This is the "signature expression reshapes over time" visual:
the cascade lights up with gene-specific TIMING (LDLR: chol->CAD->MI; BRCA: breast+cascade; TTN: HF/AF).
eta(t) IS the carrier-vs-non-carrier contrast (phi_tilde_carrier - phi_tilde_control = eta on scoped cells).
(B) TIMING ATLAS: across ALL genes, the home-disease peak-eta age (when the gene acts) -- one strip.
Contrast make_polygenic_scale_contrast.py: PRS only SCALES loading (level shift), never adds this age-structure.
INPUT : clinvar_sweep_gp/*_etat.npz (eta_t (D,T), gate (D,), inc_excess (D,T), Dpen (D,), peak_age, ages, ...)
OUTPUT: fig_monogenic_reshape_flagships.{png,pdf} + fig_monogenic_timing_atlas.{png,pdf} (+prov).
RUN via new_env_pyro2. Aggregates only (fitted eta, no EIDs).
TRIPLET COPY (2026-10-07) of make_monogenic_reshape_atlas.py (that file and ~/mendelnoulli-paper/analysis/ are left
untouched). run_reshape_atlas_triplet.sh passes SWEEPDIR / OUTSUFFIX / FLAGSHIPS / REGROW_MLH1, which the original
does NOT read -- run with them it silently read the bar sweep and overwrote the bar figures (restored 2026-10-07).
This copy reads all four: SWEEPDIR = atlas dir; OUTSUFFIX appended to every output name (required when SWEEPDIR is
set, so bar figures are never overwritten); FLAGSHIPS = comma list, PINNED (no back-fill); REGROW_MLH1 = registry
eta npz for the MLH1 row.
Panel d Registry-anchored carrier penetrance¶
Script: src/paper_keystone/Fig4_monogenic_clinvar/assemble_fig3d_triplet.py
Upstream: cancer fits: src/reparam_fig3/registry_reanchor_cancer_gp.py
Script header (inputs, recipe, outputs)
Combined Figure 6 WITH error bars.
Top row : curated non-cancer matched-control eta(t) (L0W=0, gate off) -> OBSERVED carrier KM + Greenwood 95% CI band.
Bottom row: base-free eta(t) cancer -> registry carrier Greenwood 95% CI band (same as fig_carrier_ci.png).
BOTH rows now show exactly ONE band: the 95% CI of the OBSERVED comparator (in-cohort KM top, registry KM bottom).
The model curves are drawn as lines with no band. A second, crimson binomial band used to be drawn around the
model carrier curve here (hw = 1.96*sqrt(f(1-f)/n_carr)); it was removed 2026-09-02 because (a) it is not a model
uncertainty -- it is the sampling noise n_carr people would show IF the model rate were true, so it invited the
wrong read that band-overlap = agreement, when the actual test is model-line-inside-observed-CI; (b) in the
TTN->AF panel it stayed centred on the NATIVE d['carrier'] while the thick red line is the recalibrated curve,
a figure<->code mismatch; and (c) the cancer row never had it, so the two rows were not on the same footing.
PCSK9 panel uses carriers_by_gene (same set as the OLD working panel), now under the consistent L0W=0 recipe.
Panel e Penetrance vs loading coherence¶
Script: src/clinical_app/fig_penetrance_loading_coherence_triplet.py
Upstream: age-varying atlas (run_atlas_triplet_fixed.sh)
Script header (inputs, recipe, outputs)
COHERENCE panel (Fig 5c) rebuilt on the TRIPLET backbone -- a like-for-like copy of
fig_penetrance_loading_coherence.py (the published delta0 version, left untouched).
WHAT IS HELD FIXED (identical to the published script)
genes LDLR APOB BRCA1 BRCA2 MLH1 MSH6 TTN MYBPC3 HNF4A (DSP dropped, as published)
controls 20:1, matched on sex x 5-yr age bin, np.random.RandomState(0), drawn in the same gene order
estimand Pearson r between the two curves over ages 45-80 (indices 15:50), per gene
loading theta on the EXPRESSING signature = argmax_k psi[k, home] (not the scope)
WHAT CHANGES (both curves move together -- no backbone mixing)
red inc_excess[home] from the FIXED triplet sweep clinvar_sweep_gp_<T>_fixed/<gene>_etat.npz
(healthy-row bug fixed; breast scoped to S19)
blue carrier - control mean theta, theta reconstructed batch-by-batch from the triplet's OWN trained
checkpoints (its own per-batch gamma + delta, full fit):
lambda_k = r_k(t) + zscore_batch(G37) gamma_level_k (+ t * zscore_batch(G37) gamma_slope_k) + delta_k
lambda_H = shifted disease-free curve (a=0.6, mean +1.50) + zscore_batch(G37) gamma_health + delta_H
G37 = 36 PRS residualised on 10 PCs (one global lstsq) + sex; t = arange(T)/(T-1) (the model's time_grid)
This is the reconstruction used in fig_carrier_stacks_by_backbone.py (theta_health ~0.7 at 60, as trained).
psi the triplet's own psi (mean over its 40 batches), NOT data_for_running/initial_psi_400k.pt (delta0's)
NOT AVAILABLE the red +/-1 SD ribbon: it came from a BOOTCI re-fit of the delta0 sweep; none exists for the
triplet, so the triplet panel has no ribbon (stated in the figure).
INPUTS reparam_fig3/clinvar_sweep_gp_<T>_fixed/<gene>_etat.npz ; results_slope3way/slope3way_batch_*.pt or
results_f3way40/f3way40_b*.pt ; data_for_running/{G_matrix,reference_trajectories,baselinagefamh_withpcs}
; clinvar/carriers_<gene>_clinvarPLP_wes.txt ; healthy_ref_diseasefree_smoothed.npy
OUTPUT claudefile/clinical_app/fig_penetrance_loading_coherence_<T>.{png,pdf}, penetrance_loading_coherence_<T>.csv
RUN BACKBONE=slope3way|f3way40 /opt/miniconda3/envs/new_env_pyro2/bin/python claudefile/clinical_app/fig_penetrance_loading_coherence_triplet.py
PHI group means only.