import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy import stats
import warnings
warnings.filterwarnings('ignore')

plt.rcParams.update({
    'figure.facecolor': '#0a0a14',
    'axes.facecolor': '#151525',
    'text.color': '#e0e0e0',
    'axes.labelcolor': '#e0e0e0',
    'axes.edgecolor': '#333',
    'xtick.color': '#888',
    'ytick.color': '#888',
    'figure.dpi': 100,
    'savefig.dpi': 100,
})
print('Environment ready.')

import pandas as pd

hyp_data = [
  {
    "title": "Senolytic ABT-263 Crosses BBB for Selective Neuronal Senescence Clearance",
    "composite": 0.817,
    "mech": 0.77,
    "evid": 0.59,
    "novel": 0.85,
    "feas": 0.7,
    "impact": 0.9,
    "drug": 0.56,
    "safety": 0.52,
    "comp": 0.88,
    "data": 0.63,
    "reprod": 0.59
  },
  {
    "title": "SASP-Driven Paracrine Spreading of Senescence in Hippocampal Circuits",
    "composite": 0.767,
    "mech": 0.81,
    "evid": 0.47,
    "novel": 0.9,
    "feas": 0.55,
    "impact": 0.85,
    "drug": 0.44,
    "safety": 0.41,
    "comp": 0.88,
    "data": 0.5,
    "reprod": 0.47
  },
  {
    "title": "p16INK4a-Positive Astrocyte Accumulation Drives Synaptic Loss",
    "composite": 0.733,
    "mech": 0.72,
    "evid": 0.55,
    "novel": 0.8,
    "feas": 0.65,
    "impact": 0.75,
    "drug": 0.52,
    "safety": 0.49,
    "comp": 0.78,
    "data": 0.59,
    "reprod": 0.55
  },
  {
    "title": "Microglial Senescence Impairs Amyloid-beta Phagocytosis",
    "composite": 0.783,
    "mech": 0.68,
    "evid": 0.68,
    "novel": 0.75,
    "feas": 0.8,
    "impact": 0.8,
    "drug": 0.64,
    "safety": 0.6,
    "comp": 0.78,
    "data": 0.72,
    "reprod": 0.68
  },
  {
    "title": "Dasatinib+Quercetin Combination Therapy for Brain Senescence",
    "composite": 0.783,
    "mech": 0.54,
    "evid": 0.77,
    "novel": 0.6,
    "feas": 0.9,
    "impact": 0.85,
    "drug": 0.72,
    "safety": 0.68,
    "comp": 0.72,
    "data": 0.81,
    "reprod": 0.77
  },
  {
    "title": "Senescence-Associated Mitochondrial Dysfunction in Neurons",
    "composite": 0.667,
    "mech": 0.63,
    "evid": 0.51,
    "novel": 0.7,
    "feas": 0.6,
    "impact": 0.7,
    "drug": 0.48,
    "safety": 0.45,
    "comp": 0.7,
    "data": 0.54,
    "reprod": 0.51
  },
  {
    "title": "Immune Checkpoint Modulation for Senescent Cell Clearance in CNS",
    "composite": 0.687,
    "mech": 0.79,
    "evid": 0.34,
    "novel": 0.88,
    "feas": 0.4,
    "impact": 0.78,
    "drug": 0.32,
    "safety": 0.3,
    "comp": 0.83,
    "data": 0.36,
    "reprod": 0.34
  }
]

df = pd.DataFrame(hyp_data)
df = df.rename(columns={'title': 'Hypothesis', 'composite': 'Score'})
df[['Hypothesis', 'Score', 'mech', 'evid', 'novel', 'feas', 'impact', 'drug']]

fc_data = {"CDKN2A": 2.5, "TP53": 1.8, "BCL2": 1.6, "IL6": 2.2, "CXCL1": 1.9, "MMP3": 1.4, "LMNB1": -1.7, "GLB1": 1.3}
genes = list(fc_data.keys())
np.random.seed(42)

n_samples = 20
results = []
for gene in genes:
    fc = fc_data[gene]
    control = np.random.normal(loc=8.0, scale=0.8, size=n_samples)
    disease = np.random.normal(loc=8.0 + fc, scale=1.0, size=n_samples)
    t_stat, p_val = stats.ttest_ind(control, disease)
    log2fc = np.mean(disease) - np.mean(control)
    results.append({
        'gene': gene, 'log2fc': log2fc, 'p_value': p_val,
        'neg_log10_p': -np.log10(max(p_val, 1e-10)),
        'control_mean': np.mean(control), 'disease_mean': np.mean(disease),
    })

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))

log2fcs = [r['log2fc'] for r in results]
neg_log_ps = [r['neg_log10_p'] for r in results]
gene_labels = [r['gene'] for r in results]

colors = ['#ef5350' if abs(fc) > 0.5 and nlp > 1.3 else '#888888'
          for fc, nlp in zip(log2fcs, neg_log_ps)]

ax1.scatter(log2fcs, neg_log_ps, c=colors, s=100, alpha=0.8, edgecolors='#333')
for i, gene in enumerate(gene_labels):
    ax1.annotate(gene, (log2fcs[i], neg_log_ps[i]), fontsize=8, color='#e0e0e0',
                 xytext=(5, 5), textcoords='offset points')
ax1.axhline(y=1.3, color='#ffd54f', linestyle='--', alpha=0.5, label='p=0.05')
ax1.axvline(x=-0.5, color='#888', linestyle='--', alpha=0.3)
ax1.axvline(x=0.5, color='#888', linestyle='--', alpha=0.3)
ax1.set_xlabel('log2(Fold Change)', fontsize=11)
ax1.set_ylabel('-log10(p-value)', fontsize=11)
ax1.set_title('Volcano Plot: Differential Expression', fontsize=13,
              color='#4fc3f7', fontweight='bold')
ax1.legend(fontsize=8, facecolor='#151525', edgecolor='#333', labelcolor='#e0e0e0')

x = np.arange(len(genes))
width = 0.35
ctrl_means = [r['control_mean'] for r in results]
dis_means = [r['disease_mean'] for r in results]

ax2.bar(x - width/2, ctrl_means, width, label='Control', color='#4fc3f7', alpha=0.8)
ax2.bar(x + width/2, dis_means, width, label='Disease', color='#ef5350', alpha=0.8)
ax2.set_xticks(x)
ax2.set_xticklabels(genes, rotation=45, ha='right', fontsize=9)
ax2.set_ylabel('Expression Level (log2)', fontsize=11)
ax2.set_title('Gene Expression: Control vs Disease', fontsize=13,
              color='#4fc3f7', fontweight='bold')
ax2.legend(fontsize=9, facecolor='#151525', edgecolor='#333', labelcolor='#e0e0e0')

plt.tight_layout()
plt.show()

print("\nDifferential Expression Summary")
print("=" * 70)
print(f"{'Gene':<15} {'log2FC':>10} {'p-value':>12} {'Significant':>12}")
print("-" * 70)
for r in sorted(results, key=lambda x: x['p_value']):
    sig = 'YES' if abs(r['log2fc']) > 0.5 and r['p_value'] < 0.05 else 'no'
    print(f"{r['gene']:<15} {r['log2fc']:>10.3f} {r['p_value']:>12.2e} {sig:>12}")

np.random.seed(42)
pathways = ["Cellular Senescence", "SASP Secretion", "p53/p21 Pathway", "p16INK4a/Rb Pathway", "NF-kB Inflammatory Signaling", "Senolytic Drug Targets", "Autophagy-Lysosome Pathway", "Mitochondrial Dysfunction", "DNA Damage Response", "Immune Surveillance of Senescent Cells", "BBB Permeability", "Neuroinflammation"]

enrichment_scores = np.random.exponential(2, len(pathways)) + 1
p_values = 10 ** (-np.random.uniform(1, 8, len(pathways)))
gene_counts = np.random.randint(2, 6, len(pathways))

idx = np.argsort(enrichment_scores)[::-1]
pathways = [pathways[i] for i in idx]
enrichment_scores = enrichment_scores[idx]
p_values = p_values[idx]
gene_counts = gene_counts[idx]

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 8))

sizes = gene_counts * 30
colors = -np.log10(p_values)
scatter = ax1.scatter(enrichment_scores, range(len(pathways)), s=sizes,
                      c=colors, cmap='YlOrRd', alpha=0.8, edgecolors='#333')
ax1.set_yticks(range(len(pathways)))
ax1.set_yticklabels(pathways, fontsize=9)
ax1.set_xlabel('Enrichment Score', fontsize=11)
ax1.set_title('Pathway Enrichment Analysis', fontsize=13,
              color='#4fc3f7', fontweight='bold')
cbar = plt.colorbar(scatter, ax=ax1, shrink=0.6)
cbar.set_label('-log10(p-value)', fontsize=9, color='#e0e0e0')

bar_colors = ['#ef5350' if p < 0.001 else '#ff8a65' if p < 0.01 else '#ffd54f' if p < 0.05 else '#888'
              for p in p_values]
ax2.barh(range(len(pathways)), -np.log10(p_values), color=bar_colors, alpha=0.8, edgecolor='#333')
ax2.set_yticks(range(len(pathways)))
ax2.set_yticklabels(pathways, fontsize=9)
ax2.set_xlabel('-log10(p-value)', fontsize=11)
ax2.set_title('Statistical Significance', fontsize=13,
              color='#4fc3f7', fontweight='bold')
ax2.axvline(x=-np.log10(0.05), color='#ffd54f', linestyle='--', alpha=0.7, label='p=0.05')
ax2.axvline(x=-np.log10(0.001), color='#ef5350', linestyle='--', alpha=0.7, label='p=0.001')
ax2.legend(fontsize=8, facecolor='#151525', edgecolor='#333', labelcolor='#e0e0e0')

plt.tight_layout()
plt.show()

print("\nPathway Enrichment Summary")
print("=" * 80)
print(f"{'Pathway':<40} {'Enrichment':>12} {'p-value':>12} {'Genes':>8}")
print("-" * 80)
for pw, es, pv, gc in zip(pathways, enrichment_scores, p_values, gene_counts):
    print(f"{pw:<40} {es:>12.2f} {pv:>12.2e} {gc:>8}")

hyp_data = [
  {
    "title": "Senolytic ABT-263 Crosses BBB for Selective Neuronal Senescence Clearance",
    "composite": 0.817,
    "mech": 0.77,
    "evid": 0.59,
    "novel": 0.85,
    "feas": 0.7,
    "impact": 0.9,
    "drug": 0.56,
    "safety": 0.52,
    "comp": 0.88,
    "data": 0.63,
    "reprod": 0.59
  },
  {
    "title": "SASP-Driven Paracrine Spreading of Senescence in Hippocampal Circuits",
    "composite": 0.767,
    "mech": 0.81,
    "evid": 0.47,
    "novel": 0.9,
    "feas": 0.55,
    "impact": 0.85,
    "drug": 0.44,
    "safety": 0.41,
    "comp": 0.88,
    "data": 0.5,
    "reprod": 0.47
  },
  {
    "title": "p16INK4a-Positive Astrocyte Accumulation Drives Synaptic Loss",
    "composite": 0.733,
    "mech": 0.72,
    "evid": 0.55,
    "novel": 0.8,
    "feas": 0.65,
    "impact": 0.75,
    "drug": 0.52,
    "safety": 0.49,
    "comp": 0.78,
    "data": 0.59,
    "reprod": 0.55
  },
  {
    "title": "Microglial Senescence Impairs Amyloid-beta Phagocytosis",
    "composite": 0.783,
    "mech": 0.68,
    "evid": 0.68,
    "novel": 0.75,
    "feas": 0.8,
    "impact": 0.8,
    "drug": 0.64,
    "safety": 0.6,
    "comp": 0.78,
    "data": 0.72,
    "reprod": 0.68
  },
  {
    "title": "Dasatinib+Quercetin Combination Therapy for Brain Senescence",
    "composite": 0.783,
    "mech": 0.54,
    "evid": 0.77,
    "novel": 0.6,
    "feas": 0.9,
    "impact": 0.85,
    "drug": 0.72,
    "safety": 0.68,
    "comp": 0.72,
    "data": 0.81,
    "reprod": 0.77
  },
  {
    "title": "Senescence-Associated Mitochondrial Dysfunction in Neurons",
    "composite": 0.667,
    "mech": 0.63,
    "evid": 0.51,
    "novel": 0.7,
    "feas": 0.6,
    "impact": 0.7,
    "drug": 0.48,
    "safety": 0.45,
    "comp": 0.7,
    "data": 0.54,
    "reprod": 0.51
  },
  {
    "title": "Immune Checkpoint Modulation for Senescent Cell Clearance in CNS",
    "composite": 0.687,
    "mech": 0.79,
    "evid": 0.34,
    "novel": 0.88,
    "feas": 0.4,
    "impact": 0.78,
    "drug": 0.32,
    "safety": 0.3,
    "comp": 0.83,
    "data": 0.36,
    "reprod": 0.34
  }
]

print("=" * 70)
print("STATISTICAL ANALYSIS OF HYPOTHESIS SCORES")
print("=" * 70)

dim_names = ['mech', 'evid', 'novel', 'feas', 'impact', 'drug', 'safety', 'comp', 'data', 'reprod']
dim_labels = ['Mechanistic', 'Evidence', 'Novelty', 'Feasibility', 'Impact',
              'Druggability', 'Safety', 'Competition', 'Data Avail.', 'Reproducibility']

scores_matrix = np.array([[h.get(k, 0) for k in dim_names] for h in hyp_data])

print("\n1. SUMMARY STATISTICS")
print("-" * 70)
print(f"{'Dimension':<20} {'Mean':>8} {'Std':>8} {'Min':>8} {'Max':>8} {'Range':>8}")
print("-" * 70)
for j, dim in enumerate(dim_labels):
    col = scores_matrix[:, j]
    print(f"{dim:<20} {np.mean(col):>8.3f} {np.std(col):>8.3f} "
          f"{np.min(col):>8.3f} {np.max(col):>8.3f} {np.max(col)-np.min(col):>8.3f}")

composites = [h.get('composite', 0) for h in hyp_data]
print(f"\n2. COMPOSITE SCORE DISTRIBUTION")
print("-" * 70)
print(f"Mean: {np.mean(composites):.3f}")
print(f"Median: {np.median(composites):.3f}")
print(f"Std Dev: {np.std(composites):.3f}")
if len(composites) >= 3:
    stat, p = stats.shapiro(composites)
    print(f"Shapiro-Wilk test: W={stat:.4f}, p={p:.4f} ({'Normal' if p > 0.05 else 'Non-normal'})")

if len(hyp_data) >= 4:
    top_half = scores_matrix[:len(hyp_data)//2]
    bottom_half = scores_matrix[len(hyp_data)//2:]
    print(f"\n3. TOP vs BOTTOM HYPOTHESIS COMPARISON")
    print("-" * 70)
    for j, dim in enumerate(dim_labels[:6]):
        t, p = stats.ttest_ind(top_half[:, j], bottom_half[:, j])
        sig = '*' if p < 0.05 else ''
        print(f"{dim:<20} top={np.mean(top_half[:,j]):.3f} bot={np.mean(bottom_half[:,j]):.3f} "
              f"t={t:>6.2f} p={p:.3f} {sig}")

print("\n" + "=" * 70)
print("Analysis complete.")