Reparam PRS-Signature Associations (v1.1)

Purpose. Produce the supplement-style bar plot + FDR-significant heatmap + complete heatmap (as in Extended Data Fig 9 / Fig S27) using the non-centered (reparam) γ matrices from the 400K-pooled-warm-start pool. Includes the centered (published) pool alongside for direct comparison.

Why. The published Fig S27 used the centered MAP parameterization. Referee #2 / Nature v1.1 revision asked for a sensitivity check that the headline biology holds under the non-centered (reparam) parameterization. This notebook is that check, reproducible from the CSVs deposited in docs/reparam/.

Two pools loaded (756 PRS-signature pairs each):

Pool	CSV	γ definition	FDR<0.05 cells
Centered (published)	prs_signature_centered_gamma_associations.csv	λ ~ GP(μ + Gγ, K), λ free	116
Reparam 400K-init	prs_signature_400k_gamma_associations.csv	λ = μ + Gγ + δ, δ ~ GP(0, K)	416

Both CSVs were generated by generate_prs_signature_plots.py (Z = γ̄/SEM across 40 leave-one-out batches; Benjamini–Hochberg FDR across 756 pairs).

Headline result. The named biological hits (CAD→Sig5, T2D→Sig15, LDL→Sig5, BMI→Sig15, AF→Sig0, HT→Sig5) have the same sign and strong significance in both pools.

import sys
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from IPython.display import display
from matplotlib.patches import Patch

HERE = Path('.').resolve()
if HERE.name != 'reparam':
    # Notebook may be opened from repo root; resolve relative to known docs/reparam path.
    HERE = Path('/Users/sarahurbut/aladynoulli2/docs/reparam')

C = pd.read_csv(HERE / 'prs_signature_centered_gamma_associations.csv')
N = pd.read_csv(HERE / 'prs_signature_400k_gamma_associations.csv')

print(f'Centered (published) rows: {len(C)}, FDR<0.05: {C.significant_fdr.sum()}')
print(f'Reparam 400K-init    rows: {len(N)}, FDR<0.05: {N.significant_fdr.sum()}')

# Disease-category colours (same scheme as published Fig S27)
CATEGORY_COLORS = {
    'Cardiovascular': '#d62728',
    'Metabolic':      '#1f77b4',
    'Neurological':   '#2ca02c',
    'Cancer':         '#9467bd',
    'Autoimmune':     '#ff7f0e',
    'Other':          '#7f7f7f',
}

Centered (published) rows: 756, FDR<0.05: 116
Reparam 400K-init    rows: 756, FDR<0.05: 416

Plot helpers (mirror generate_prs_signature_plots.py)

Three plots per pool: (1) top-30 PRS-signature pairs by |Z|, coloured by disease category; (2) FDR<0.05 heatmap; (3) complete 36×21 Z-score heatmap.

def plot_top_bar(df, label, n_top=30):
    df = df.copy()
    df['absz'] = df['z_score'].abs()
    top = df.dropna(subset=['z_score']).nlargest(n_top, 'absz').sort_values('z_score')
    fig, ax = plt.subplots(figsize=(10, max(8, len(top) * 0.35)))
    colors = [CATEGORY_COLORS.get(c, CATEGORY_COLORS['Other']) for c in top['category']]
    ax.barh(range(len(top)), top['z_score'], color=colors, alpha=0.85)
    ax.set_yticks(range(len(top)))
    def make_label(r):
        s = f"{r['prs']} - {r['signature']}"
        if not np.isnan(r['p_value_fdr']):
            if r['p_value_fdr'] < 0.001:
                s += ' (FDR<0.001)'
            else:
                s += f" (FDR={r['p_value_fdr']:.3f})"
            if r['significant_fdr']:
                s += '*'
        return s
    ax.set_yticklabels(top.apply(make_label, axis=1), fontsize=9)
    ax.set_xlabel('Z-score (effect / SEM)', fontsize=12, fontweight='bold')
    ax.set_title(f'Top {n_top} PRS-Signature Associations -- {label}', fontsize=13, fontweight='bold', pad=10)
    ax.axvline(x=0, color='black', linestyle='--', linewidth=0.5, alpha=0.5)
    ax.grid(True, alpha=0.3, axis='x')
    handles = [Patch(facecolor=c, label=cat, alpha=0.85) for cat, c in CATEGORY_COLORS.items()]
    ax.legend(handles=handles, loc='lower right', fontsize=9, framealpha=0.95)
    plt.tight_layout()
    plt.show()


def plot_heatmap(df, label, fdr_only=False):
    sub = df[df['p_value_fdr'] < 0.05].copy() if fdr_only else df.copy()
    pivot = sub.pivot_table(values='z_score', index='prs', columns='signature', aggfunc='mean')
    # Order PRS by category then alphabetically
    cat_map = dict(zip(df['prs'].astype(str), df['category']))
    order = []
    for cat in ['Cardiovascular', 'Metabolic', 'Autoimmune', 'Neurological', 'Cancer', 'Other']:
        order.extend(sorted([p for p in pivot.index if cat_map.get(str(p)) == cat]))
    pivot = pivot.reindex([p for p in order if p in pivot.index])
    fig, ax = plt.subplots(figsize=(max(12, len(pivot.columns)*0.8), max(10, len(pivot)*0.4)))
    sns.heatmap(pivot, cmap='RdBu_r', center=0, vmin=-5, vmax=5,
                cbar_kws={'label': 'Z-score'}, ax=ax, linewidths=0.5, linecolor='gray')
    suffix = ' (FDR<0.05)' if fdr_only else ' (all 36 PRS x 21 signatures)'
    ax.set_title(f'PRS-Signature Z-scores -- {label}{suffix}', fontsize=13, fontweight='bold', pad=10)
    ax.set_xlabel('Signature', fontsize=12, fontweight='bold')
    ax.set_ylabel('Polygenic Risk Score', fontsize=12, fontweight='bold')
    plt.xticks(rotation=45, ha='right')
    plt.yticks(rotation=0)
    plt.tight_layout()
    plt.show()

Centered (published) — three plots

These are equivalent to the published Fig S27 panels A, B, C, drawn directly from the deposited CSV.

plot_top_bar(C, 'Centered (published)')
plot_heatmap(C, 'Centered (published)', fdr_only=True)
plot_heatmap(C, 'Centered (published)', fdr_only=False)

Reparam 400K-init (v1.1) — three plots

Same plot code applied to the reparam-400K γ matrix. Both pools are aggregated across 40 LOO batches; reparam pool warm-starts γ from a shared 400K-patient signature-score ridge regression.

plot_top_bar(N, 'Reparam 400K-init (v1.1)')
plot_heatmap(N, 'Reparam 400K-init (v1.1)', fdr_only=True)
plot_heatmap(N, 'Reparam 400K-init (v1.1)', fdr_only=False)

Headline named hits — side-by-side

The biological hits the paper text highlights, in both pools.

named = [('CAD', 'Sig 5'), ('LDL_SF', 'Sig 5'), ('T2D', 'Sig 15'),
         ('BMI', 'Sig 15'), ('AF', 'Sig 0'), ('HT', 'Sig 5'), ('CVD', 'Sig 5')]
m = (C.rename(columns=lambda c: c+'_c' if c not in ('prs','signature') else c)
     .merge(N.rename(columns=lambda c: c+'_n' if c not in ('prs','signature') else c),
            on=['prs','signature']))
rows = []
for p, s in named:
    r = m[(m['prs']==p) & (m['signature']==s)]
    if not len(r): continue
    r = r.iloc[0]
    rows.append({
        'PRS': p, 'Signature': s,
        'centered_gamma': f'{r.effect_c:+.4f}',  'centered_Z': f'{r.z_score_c:+.1f}',
        '400k_gamma':     f'{r.effect_n:+.4f}',  '400k_Z':     f'{r.z_score_n:+.1f}',
        'sign_agree': int(np.sign(r.effect_c) == np.sign(r.effect_n)),
    })
display(pd.DataFrame(rows))

	PRS	Signature	centered_gamma	centered_Z	400k_gamma	400k_Z	sign_agree
0	CAD	Sig 5	+0.1534	+27.2	+0.3345	+35.8	1
1	LDL_SF	Sig 5	+0.0713	+22.7	+0.1380	+13.7	1
2	T2D	Sig 15	+0.1539	+58.3	+1.1458	+79.1	1
3	BMI	Sig 15	+0.0220	+13.7	+0.2309	+18.0	1
4	AF	Sig 0	+0.0379	+22.1	+0.3006	+35.1	1
5	HT	Sig 5	+0.0475	+12.3	+0.2088	+26.4	1
6	CVD	Sig 5	+0.0770	+14.2	+0.1687	+22.0	1

Z-score scatter (consistency check)

Correlation of the 756-cell Z map between the two pools, restricted to cells where both have a defined Z (Sig 20 is the reference and has no γ → 36 cells dropped).

mask = m.z_score_c.notna() & m.z_score_n.notna()
x = m.loc[mask, 'z_score_c'].values
y = m.loc[mask, 'z_score_n'].values
fig, ax = plt.subplots(figsize=(7, 7))
ax.scatter(x, y, s=12, alpha=0.55, c='#1a3a5c')
lim = max(np.abs(x).max(), np.abs(y).max()) * 1.05
ax.plot([-lim, lim], [-lim, lim], '--', color='red', alpha=0.5, lw=1)
ax.axhline(0, color='gray', lw=0.5); ax.axvline(0, color='gray', lw=0.5)
ax.set_xlabel('Z (Centered / published)')
ax.set_ylabel('Z (Reparam 400K-init / v1.1)')
r = np.corrcoef(x, y)[0, 1]
strong = (np.abs(x) > 2) | (np.abs(y) > 2)
r_s = np.corrcoef(x[strong], y[strong])[0, 1] if strong.sum() > 1 else np.nan
ax.set_title(f'r(all) = {r:.3f}   r(|Z|>2) = {r_s:.3f}    (red dashed = y=x)', fontsize=11)
ax.set_xlim(-lim, lim); ax.set_ylim(-lim, lim)
plt.tight_layout()
plt.show()

Provenance

• Centered γ source: paper_figs/fig4/gamma_adventures/prs_signatures_nolr/gamma_associations.csv ≡ SuppDataFiles-3/Annotation/gamma_associations.csv (bit-for-bit). 40 LOO batches of UKB centered MAP fit; γ enters via GP prior mean of free λ.
• Reparam 400K-init γ source: 40 LOO batches at Dropbox/censor_e_batchrun_vectorized_REPARAM_v3_nokappa_400k_init/. Each batch warm-starts γ from gamma_level_pooled_400k.pt (signature-score ridge regression on 400K UKB participants). Trainer: claudefile/train_nokappa_v3_all40_400k_init.py.
• Z / FDR computation: pyScripts/dec_6_revision/new_notebooks/main_paper_figures/generate_prs_signature_plots.py (Z = γ̄ / SEM across 40 batches; BH-FDR across all 756 pairs).
• CSVs in this folder: prs_signature_centered_gamma_associations.csv, prs_signature_400k_gamma_associations.csv.
• Comparison driver script: claudefile/compare_prs_signature_pools.py.

Conclusion

Both pools agree on the direction and strong significance of every headline biological hit. The reparam pool finds ~3.6× more FDR-significant cells overall (γ has full likelihood gradient and isn't damped by the GP prior), but those additional cells are weak signals near the noise floor — not new biology. No revision to Fig S27 is required; the reparam pool is a clean v1.1 sensitivity check.