"""
Supplementary Figure 35 — MCF-7 per-chromosome periodicity strength: chromosome-level
coarse features (peak abundance, dominant period) do not predict sig95.

REAL DATA ONLY. No rng. Replaces the v21 synthetic placeholder
(`extended_data_figure_9_cnv.py`, rng.uniform(0.25, 0.55, n_chroms)) which fabricated
sig95 values per chromosome.

Source CSV (BIOMNI 2026-04-04 ENCODE DNase pan-cell-line sweep):
  - encode_dnase_wii_per_chromosome.csv (10 cell lines × 22 autosomes; manuscript-canonical
    chromperiod Paul m=2 pipeline applied to ENCODE DNase-seq narrowPeak)

The figure addresses the methodological question "could the wave just reflect peak-density
or per-chromosome peak abundance?" using real MCF-7 data alone:

  - Panel a: Per-chromosome sig95 across 22 MCF-7 autosomes. chr17 (the lowest, sig95 = 18.9%)
    is flagged because chr17q carries the canonical HER2 amplicon in MCF-7
    (literature-attested), and chr20 (second-lowest at 23.3%) carries the 20q13 amplicon.
    The point is honest disclosure of two known MCF-7 structural-rearrangement hotspots
    against the genome-wide range 19-36%.
  - Panel b: sig95 vs peak count per chromosome. Spearman ρ tests whether periodicity strength
    is just a peak-density artefact. Real data: sig95 is essentially uncorrelated with peak
    count — the wave is not a peak-density readout.
  - Panel c: sig95 vs dominant period (Mbp). Shows that periodicity strength and wavelength
    are internally consistent (longer chromosomes have longer dominant periods, but sig95
    spans 19-36% across both short and long chromosomes).

Caption claim: "Coarse chromosome-level features do not predict periodicity strength in
MCF-7"; the chr17/chr20 outliers are honestly flagged as known amplicon territories.
"""
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-relative paths -----
BUNDLE_ROOT = Path(__file__).resolve().parents[2]
DATA_DIR = BUNDLE_ROOT / '04_supplementary_data' / 'figure_data'
SRC_CSV = DATA_DIR / 'encode_dnase_wii_per_chromosome.csv'

OUT_BASENAMES = [
    BUNDLE_ROOT / '02_figures' / 'supplementary' / 'supplementary_figure_35_mcf7_periodicity',
    BUNDLE_ROOT / '02_figures' / 'extended_data' / 'extended_data_figure_9_cnv',
    BUNDLE_ROOT / '03_figures_to_edit' / 'supplementary' / 'supplementary_figure_35_mcf7_periodicity',
    BUNDLE_ROOT / '03_figures_to_edit' / 'extended_data' / 'extended_data_figure_9_cnv',
]
DRAFT_BASENAME = BUNDLE_ROOT / '05_code' / 'figure_regen' / '_drafts' / 'supp35_mcf7_v1'

# ----- Load real data -----
df = pd.read_csv(SRC_CSV)
mcf = df[df['cell_line'] == 'MCF-7'].copy()
# Sort chromosomes naturally
mcf['_k'] = mcf['chrom'].str.replace('chr', '').astype(int)
mcf = mcf.sort_values('_k').drop(columns='_k').reset_index(drop=True)

print(f"[verify] MCF-7 autosomes loaded: {len(mcf)}")
print(f"[verify] sig95 range: {mcf['sig95'].min():.2f}% – {mcf['sig95'].max():.2f}%")
print(f"[verify] median sig95 = {mcf['sig95'].median():.2f}%")
print(f"[verify] chr17 sig95 = {mcf.loc[mcf['chrom']=='chr17','sig95'].values[0]:.2f}% (HER2-amplicon-bearing)")
print(f"[verify] chr20 sig95 = {mcf.loc[mcf['chrom']=='chr20','sig95'].values[0]:.2f}% (20q13-amplicon)")

# Spearman: sig95 vs peak count
rho_pk, p_pk = spearmanr(mcf['n_peaks'], mcf['sig95'])
print(f"[verify] sig95 vs n_peaks: Spearman ρ = {rho_pk:+.3f}, p = {p_pk:.3f}")

# Spearman: sig95 vs dominant period
rho_T, p_T = spearmanr(mcf['dom_period_Mbp'], mcf['sig95'])
print(f"[verify] sig95 vs dom_period_Mbp: Spearman ρ = {rho_T:+.3f}, p = {p_T:.3f}")

# Mann-Whitney: chr17+chr20 (amplified) vs the rest
from scipy.stats import mannwhitneyu
amplified_mask = mcf['chrom'].isin(['chr17', 'chr20'])
amp_vals = mcf.loc[amplified_mask, 'sig95'].values
oth_vals = mcf.loc[~amplified_mask, 'sig95'].values
mw_stat, mw_p = mannwhitneyu(amp_vals, oth_vals, alternative='less')
print(f"[verify] amplified (chr17, chr20) vs others Mann-Whitney p (one-sided less) = {mw_p:.3f}")

# ============== FIGURE ==============
apply_style()
fig = plt.figure(figsize=(13.0, 4.4))
gs = fig.add_gridspec(1, 3, wspace=0.34, left=0.06, right=0.985, top=0.84, bottom=0.18,
                      width_ratios=[1.4, 1.0, 1.0])
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: per-chromosome sig95 with amplicon highlights -----
chrom_idx = np.arange(len(mcf))
chrom_labels = [c.replace('chr', '') for c in mcf['chrom'].values]
sig95 = mcf['sig95'].values
amplified = mcf['chrom'].isin(['chr17', 'chr20']).values

# Bar colours
bar_colors = [WONG['vermil'] if amplified[i] else WONG['blue'] for i in chrom_idx]

bars = ax_a.bar(chrom_idx, sig95, color=bar_colors,
                edgecolor='white', linewidth=0.6, alpha=0.92)

# Median reference line
median_sig95 = float(np.median(sig95))
ax_a.axhline(median_sig95, color=WONG['gold'], lw=1.0, ls='--', alpha=0.85,
             label=f'median = {median_sig95:.1f}%')
# AR(1) 5% null reference
ax_a.axhline(5, color=WONG['grey'], lw=0.7, ls=':', alpha=0.7,
             label='AR(1) 5% null')

# Annotate the two amplified chromosomes (place text outside bar region, above)
y_top = sig95.max() * 1.18
for c, label, dx in [('chr17', 'chr17 (17q12 / HER2)', -0.5),
                     ('chr20', 'chr20 (20q13)', +0.5)]:
    i = list(mcf['chrom'].values).index(c)
    ax_a.annotate(label,
                  xy=(i, sig95[i]),
                  xytext=(i + dx, sig95.max() * 1.10),
                  fontsize=7.5, ha=('right' if dx < 0 else 'left'), va='bottom',
                  color=WONG['vermil'], fontweight='bold',
                  arrowprops=dict(arrowstyle='->', color=WONG['vermil'], lw=0.7,
                                  shrinkA=0, shrinkB=2,
                                  connectionstyle=f"arc3,rad={0.18 * np.sign(dx)}"))

ax_a.set_xticks(chrom_idx)
ax_a.set_xticklabels(chrom_labels, fontsize=7)
ax_a.set_xlabel('MCF-7 autosome', fontsize=9)
ax_a.set_ylabel('CWT periodicity strength sig95 (%)', fontsize=9)
ax_a.set_ylim(0, sig95.max() * 1.25)
ax_a.legend(loc='upper left', fontsize=8, framealpha=0.92)
ax_a.set_title(f'All 22 autosomes above AR(1) 5% null; canonical MCF-7 amplicons sit lowest',
               loc='left', pad=8, fontsize=10)
panel_letter(fig, ax_a, 'a')

# ----- Panel b: sig95 vs peak count -----
ax_b.scatter(mcf['n_peaks'].values, sig95,
             s=46, c=bar_colors, edgecolor='white', linewidth=0.5, alpha=0.92, zorder=3)
# Annotate amplified chromosomes
for i, c in enumerate(mcf['chrom'].values):
    if amplified[i]:
        ax_b.annotate(c.replace('chr', 'chr'),
                      xy=(mcf['n_peaks'].values[i], sig95[i]),
                      xytext=(4, 4), textcoords='offset points',
                      fontsize=7, color=WONG['vermil'], fontweight='bold')
# Trend line (visual only)
zfit = np.polyfit(mcf['n_peaks'].values, sig95, 1)
xs = np.array([mcf['n_peaks'].min(), mcf['n_peaks'].max()])
ax_b.plot(xs, np.polyval(zfit, xs), color=WONG['grey'], lw=0.8, ls='--', alpha=0.7)
ax_b.text(0.04, 0.96,
          f'Spearman ρ = {rho_pk:+.2f}\np = {p_pk:.2f}\nn = 22 autosomes',
          transform=ax_b.transAxes, fontsize=8, va='top', ha='left',
          bbox=dict(boxstyle='round,pad=0.3', facecolor='white',
                    edgecolor='0.7', lw=0.4, alpha=0.95))
ax_b.set_xlabel('peaks per chromosome (DNase narrowPeak count)', fontsize=9)
ax_b.set_ylabel('sig95 (%)', fontsize=9)
ax_b.set_title('Amplicon-bearing chromosomes sit below the trend',
               loc='left', pad=8, fontsize=10)
panel_letter(fig, ax_b, 'b')

# ----- Panel c: sig95 vs dominant period -----
ax_c.scatter(mcf['dom_period_Mbp'].values, sig95,
             s=46, c=bar_colors, edgecolor='white', linewidth=0.5, alpha=0.92, zorder=3)
# Annotate amplified chromosomes
for i, c in enumerate(mcf['chrom'].values):
    if amplified[i]:
        ax_c.annotate(c.replace('chr', 'chr'),
                      xy=(mcf['dom_period_Mbp'].values[i], sig95[i]),
                      xytext=(4, 4), textcoords='offset points',
                      fontsize=7, color=WONG['vermil'], fontweight='bold')
ax_c.text(0.04, 0.96,
          f'Spearman ρ = {rho_T:+.2f}\np = {p_T:.2f}\nn = 22',
          transform=ax_c.transAxes, fontsize=8, va='top', ha='left',
          bbox=dict(boxstyle='round,pad=0.3', facecolor='white',
                    edgecolor='0.7', lw=0.4, alpha=0.95))
# amplified-vs-others MW p
ax_c.text(0.04, 0.74,
          f'amplified < others:\nMW $p = {mw_p:.3f}$',
          transform=ax_c.transAxes, fontsize=8, va='top', ha='left',
          color=WONG['vermil'], fontweight='bold',
          bbox=dict(boxstyle='round,pad=0.3', facecolor='white',
                    edgecolor=WONG['vermil'], lw=0.5, alpha=0.95))
ax_c.set_xlabel('dominant period $T$ (Mbp)', fontsize=9)
ax_c.set_ylabel('sig95 (%)', fontsize=9)
ax_c.set_title(f'Wavelength × strength: chr17 + chr20 fall off the trend',
               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] {base}.png + .svg")
else:
    print("[draft only — set PROMOTE=1 to write canonical paths]")
plt.close(fig)
