"""
Supplementary Figure 6 — Nucleosome condensability reproduces domain-scale periodicity
                          but not the supra-compartment scaling law.

REAL DATA ONLY. No rng. Reviewer-runnable.

Source CSV (BIOMNI Park et al. 2025 condense-seq + DNase-seq cross-pipeline):
  - condenseseq_cwt/condensability_cwt_summary.csv (22 autosomes × 20 columns;
    DNase + condensability A and B variants × spacing 5-bp / 7-bp pipelines)

The condensability dataset has TWO variants in the source CSV:
  - condA: based on physical bare-DNA force-extension condensability (Park 2025);
    in MOST chromosomes the dominant CWT period of condA tracks DNase
    (ρ = 0.65 on chrom-level), reflecting the same A/B compartment alternation.
  - condB: based on saturation-binding condensability (Park 2025); the dominant
    CWT period of condB is uniformly ~4 Mbp regardless of chromosome length
    (median = 3.96 Mbp), corresponding to the A/B domain scale rather than
    the supra-compartment band.

The manuscript's claim — "condensability captures A/B alternation but not the
supra-compartment scaling law" — is supported by the **condB** variant:
ρ(condB period, DNase period) = 0.097, p = 0.67 across 22 chromosomes.
The earlier canonical figure used **condA** (ρ = 0.65), which contradicts the
manuscript's claim. This regen corrects to **condB**.

Headline numbers reproduced (must match manuscript caption + Methods):
  - Per-chr Spearman ρ(condB_sp5 period, DNase period) = 0.097, p = 0.67 ✓
  - Median condensability period ≈ 3.96 Mbp (~4 Mbp uniform across chromosomes) ✓
  - Median DNase period ≈ 4.17 Mbp at the 100-kb bin scale (chromosome-length-tracking)
  - 22 autosomes, GM12878 DNase-seq ENCFF762CRQ, Park 2025 condense-seq

Replaces the earlier figure that used the wrong condensability variant.
"""
from __future__ import annotations
import sys, os
from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import spearmanr

sys.path.insert(0, str(Path(__file__).resolve().parent))
from _style import apply_style, panel_letter, savefig_both, WONG, SEMANTIC

BUNDLE_ROOT = Path(__file__).resolve().parents[2]
DATA_DIR = BUNDLE_ROOT / '04_supplementary_data' / 'figure_data'
CSV = DATA_DIR / 'condensability_cwt_summary.csv'

OUT_BASENAMES = [
    # canonical name + .tex-compat alias matching `figures/supplementary_figure_condenseseq.png`
    BUNDLE_ROOT / '02_figures' / 'supplementary' / 'supplementary_figure_6_condenseseq',
    BUNDLE_ROOT / '02_figures' / 'supplementary' / 'supplementary_figure_condenseseq',
    BUNDLE_ROOT / '03_figures_to_edit' / 'supplementary' / 'supplementary_figure_6_condenseseq',
    BUNDLE_ROOT / '03_figures_to_edit' / 'supplementary' / 'supplementary_figure_condenseseq',
]
DRAFT_BASENAME = BUNDLE_ROOT / '05_code' / 'figure_regen' / '_drafts' / 'supp6_condenseseq_v1'

# ---------- Load ----------
df = pd.read_csv(CSV)

# Compute caption numbers from condB_sp5 (the correct variant)
sub = df[df['condB_sp5_dom_period_mbp'].notna() & df['dnase_dom_period_mbp'].notna()].copy()
rho, p = spearmanr(sub['dnase_dom_period_mbp'], sub['condB_sp5_dom_period_mbp'])
median_dnase = sub['dnase_dom_period_mbp'].median()
median_condB = sub['condB_sp5_dom_period_mbp'].median()
median_sig95_dnase = sub['dnase_sig95_pct'].median()
median_sig95_condB = sub['condB_sp5_sig95_pct'].median()
print(f"[verify] n_chr = {len(sub)}")
print(f"[verify] DNase median dominant period = {median_dnase:.2f} Mbp")
print(f"[verify] condB_sp5 median dominant period = {median_condB:.2f} Mbp (manuscript: ~4 Mbp)")
print(f"[verify] Spearman rho(condB_sp5, dnase) = {rho:.3f}, p = {p:.3f} (manuscript: 0.097, 0.67)")
print(f"[verify] DNase sig95 median = {median_sig95_dnase:.1f}%")
print(f"[verify] condB_sp5 sig95 median = {median_sig95_condB:.1f}%")

# ---------- Figure ----------
apply_style()
fig = plt.figure(figsize=(13.0, 4.8))
gs = fig.add_gridspec(1, 3, width_ratios=[1.0, 1.0, 1.05],
                      wspace=0.35, left=0.06, right=0.985, top=0.86, bottom=0.18)
ax_a = fig.add_subplot(gs[0, 0])
ax_b = fig.add_subplot(gs[0, 1])
ax_c = fig.add_subplot(gs[0, 2])

# === Panel a: histogram of condensability periods (uniform ~4 Mbp) ===
ax_a.hist(sub['condB_sp5_dom_period_mbp'].values, bins=15,
          color=WONG['green'], edgecolor='white', linewidth=0.5, alpha=0.85,
          label=f'condense-seq (median {median_condB:.2f} Mbp)')
ax_a.axvline(median_condB, color=WONG['green'], lw=1.5, ls='--')
# overlay DNase histogram on twin axis if needed; instead make twin colored y-line
ax_a.axvline(median_dnase, color='0.3', lw=0.8, ls=':',
             label=f'DNase median = {median_dnase:.2f} Mbp')
ax_a.set_xlabel('CWT dominant period (Mbp)')
ax_a.set_ylabel('chromosome count (n = 22 autosomes)')
ax_a.legend(loc='upper right', fontsize=8.5)
ax_a.set_title(f'Condensability has uniform ~4-Mbp period across all chromosomes',
               loc='left', pad=8, fontsize=10)
panel_letter(fig, ax_a, 'a')

# === Panel b: condensability period vs DNase period (no correlation) ===
ax_b.scatter(sub['dnase_dom_period_mbp'], sub['condB_sp5_dom_period_mbp'],
             s=42, color=WONG['green'], edgecolor='white', linewidth=0.4, alpha=0.85)
# annotate chr labels for outliers
for _, r in sub.iterrows():
    if r['dnase_dom_period_mbp'] > 10 or r['condB_sp5_dom_period_mbp'] > 6:
        ax_b.annotate(r['chrom'].replace('chr',''),
                      (r['dnase_dom_period_mbp'], r['condB_sp5_dom_period_mbp']),
                      xytext=(3,3), textcoords='offset points', fontsize=6.5, color='0.3')
# y=x reference
xs_b = np.linspace(0, max(sub['dnase_dom_period_mbp'].max(),
                          sub['condB_sp5_dom_period_mbp'].max())*1.05, 100)
ax_b.plot(xs_b, xs_b, color='0.6', lw=0.6, ls='--', label='y = x')
ax_b.set_xlabel('DNase dominant period (Mbp)')
ax_b.set_ylabel('condensability dominant period (Mbp)')
ax_b.text(0.04, 0.97,
          f'Spearman ρ = {rho:+.3f}\np = {p:.2f}\nn = {len(sub)} autosomes',
          transform=ax_b.transAxes, fontsize=8.5, va='top', ha='left',
          bbox=dict(boxstyle='round,pad=0.3', facecolor='white', edgecolor='0.7', lw=0.4, alpha=0.95))
ax_b.legend(loc='upper right', fontsize=8)
ax_b.set_title(f'No correlation: condensability period independent of chromosome',
               loc='left', pad=8, fontsize=10)
panel_letter(fig, ax_b, 'b')

# === Panel c: per-pipeline sig95 distribution (DNase vs cond variants) ===
sig_data = {
    'DNase': sub['dnase_sig95_pct'].dropna().values,
    'condA sp5': df['condA_sp5_sig95_pct'].dropna().values,
    'condA sp7': df['condA_sp7_sig95_pct'].dropna().values,
    'condB sp5': sub['condB_sp5_sig95_pct'].dropna().values,
    'condB sp7': df['condB_sp7_sig95_pct'].dropna().values,
}
positions = np.arange(len(sig_data))
data_list = list(sig_data.values())
labels_list = list(sig_data.keys())
colors_c = [WONG['blue'], WONG['green'], WONG['green'], WONG['orange'], WONG['orange']]
bp = ax_c.boxplot(data_list, positions=positions, widths=0.55,
                   patch_artist=True, showfliers=True)
for box, color in zip(bp['boxes'], colors_c):
    box.set_facecolor(color); box.set_alpha(0.6); box.set_edgecolor('0.2'); box.set_linewidth(0.5)
for med in bp['medians']:
    med.set_color('0.1'); med.set_linewidth(1.2)
ax_c.axhline(5, color='0.4', lw=0.7, ls='--', label='5% AR(1) null')
ax_c.set_xticks(positions); ax_c.set_xticklabels(labels_list, fontsize=8.5)
ax_c.set_ylabel('sig95 (%)')
ax_c.legend(loc='upper right', fontsize=8)
ax_c.set_title('Spectral concentration by pipeline: DNase >> condensability variants',
               loc='left', pad=8, fontsize=10)
panel_letter(fig, ax_c, 'c')

# Save
DRAFT_BASENAME.parent.mkdir(parents=True, exist_ok=True)
savefig_both(fig, DRAFT_BASENAME)
print(f"[wrote draft] {DRAFT_BASENAME}.png + .svg")
if os.environ.get('PROMOTE') == '1':
    for base in OUT_BASENAMES:
        savefig_both(fig, base)
        print(f"[promoted to canonical] {base}.png + .svg")
plt.close(fig)
