はじめに

がんの診断・予後予測において、遺伝子発現データを用いた機械学習分類は臨床応用が進んでいる重要な分野です。しかし、探索的データ解析（EDA）、特徴量選択、次元削減、機械学習分類、生存解析と多段階にわたるパイプラインを一から構築するのは、バイオインフォマティクス初学者にとってハードルが高いのが現状です。

本記事では、Claude Code に自然言語で指示するだけで、scikit-learn の乳がんデータセット（569 サンプル、30 特徴量）を使ったがん遺伝子発現データの包括的分類パイプラインと、合成データによる生存解析を構築・実行した結果を報告します。

何ができるようになるか

この記事を通じて、以下のエンドツーエンドのがん分類・生存解析パイプラインを実行できるようになります。

環境

項目	バージョン
OS	Ubuntu (WSL2)
Python	3.12.3
Claude Code	最新版 (Claude Opus 4.6)
scikit-learn	1.8.0
XGBoost	3.2.0
matplotlib	3.10.8
seaborn	0.13.2
scipy	1.16.3
umap-learn	0.5.11
lifelines	0.30.1

使用データ

項目	内容
データセット	scikit-learn Breast Cancer Wisconsin (Diagnostic)
説明	FNA（穿刺吸引細胞診）画像から計算された 30 の数値特徴量
サンプル数	569 サンプル（悪性 212、良性 357）
特徴量数	30 features（mean / SE / worst の各 10 特徴量）
アクセス	`sklearn.datasets.load_breast_cancer()` で直接ロード
生存解析	合成データ（200 患者、ER+/ER-/HER2+ サブタイプ）

Step 1: 環境構築とデータロード

1-1. パッケージのインストール

# プロジェクトディレクトリの作成
mkdir -p ~/ScientificSkills/Exp-03/{data,figures,results}
cd ~/ScientificSkills/Exp-03

# 必要パッケージのインストール
pip install --break-system-packages scikit-learn matplotlib seaborn scipy umap-learn lifelines adjustText xgboost

1-2. データの読み込みと探索的データ解析（EDA）

from sklearn.datasets import load_breast_cancer
import pandas as pd

data = load_breast_cancer()
X = pd.DataFrame(data.data, columns=data.feature_names)
y = pd.Series(data.target, name="diagnosis")
target_names = data.target_names  # ['malignant', 'benign']

print(f"Dataset shape: {X.shape}")
print(f"Classes: malignant={sum(y==0)}, benign={sum(y==1)}")

Dataset shape: (569, 30)
Classes: malignant=212, benign=357

sklearn の乳がんデータセットでは 0 = 悪性（malignant）、1 = 良性（benign） として定義されています。30 個の特徴量は、細胞核画像から計算された radius、texture、perimeter、area、smoothness、compactness、concavity、concave points、symmetry、fractal dimension の mean / SE / worst の 3 カテゴリで構成されています。

1-3. 特徴量分布の可視化

悪性（赤）と良性（青）で特徴量の分布を比較すると、mean radius、mean perimeter、mean area において明確な分離が見られます。悪性腫瘍は良性に比べてこれらの値が大きく、細胞核のサイズが診断の重要な指標であることが示唆されます。

1-4. 相関ヒートマップ

分散の大きい上位 15 特徴量間の Pearson 相関を可視化しました。worst area と worst perimeter（r = 0.99）、mean area と mean perimeter（r = 0.99）など、幾何学的に関連する特徴量間に非常に強い相関が見られます。これは多重共線性の存在を示しており、正則化やPCA による次元削減が有効であることを示唆しています。

Step 2: 特徴量選択と差次発現解析

2-1. Mann-Whitney U 検定による特徴量の統計的評価

悪性 vs 良性の各特徴量について、ノンパラメトリックな Mann-Whitney U 検定 を実施し、log2 fold change と合わせてボルケーノプロットを作成しました。

from scipy.stats import mannwhitneyu

results = []
for feat in X.columns:
    mal_vals = X.loc[y == 0, feat]
    ben_vals = X.loc[y == 1, feat]
    stat, pval = mannwhitneyu(mal_vals, ben_vals, alternative="two-sided")
    fc = mal_vals.mean() / max(ben_vals.mean(), 1e-10)
    log2fc = np.log2(fc) if fc > 0 else 0
    results.append({"feature": feat, "log2_fold_change": log2fc, "p_value": pval})

de_df = pd.DataFrame(results).sort_values("p_value")
de_df["p_adjusted"] = np.minimum(de_df["p_value"] * len(de_df), 1.0)  # Bonferroni

上位 10 特徴量

特徴量	log2FC	p-value	p-adjusted
worst perimeter	0.700	2.58e-80	7.75e-79
worst radius	0.660	1.14e-78	3.41e-77
worst area	1.348	1.80e-78	5.41e-77
worst concave points	1.292	1.86e-77	5.59e-76
mean concave points	1.775	1.01e-76	3.02e-75
mean perimeter	0.563	3.55e-71	1.07e-69
mean area	1.080	1.54e-68	4.62e-67
mean concavity	1.804	2.16e-68	6.49e-67
mean radius	0.524	2.69e-68	8.08e-67
area error	1.782	5.77e-65	1.73e-63

30 特徴量すべてが Bonferroni 補正後でも p < 0.05 であり、このデータセットの判別力の高さがわかります。特に concave points（凹点の数）と concavity（凹度）は log2FC が大きく、悪性腫瘍の細胞核の形態的特徴を強く反映しています。

2-2. ボルケーノプロット

赤点は悪性で上昇（log2FC > 0.5, padj < 0.05）、青点は悪性で低下する特徴量を示します。全体的に右側（悪性で上昇）に偏っており、悪性腫瘍の細胞核は サイズ（area, perimeter, radius） と 形状の複雑さ（concavity, concave points） の両方で良性より高い値を示すことが確認できます。

Step 3: 次元削減（PCA、UMAP、t-SNE）

3-1. 3 手法の比較

標準化後のデータに対して、PCA、UMAP、t-SNE の 3 手法で 2 次元に射影しました。

from sklearn.decomposition import PCA
from sklearn.manifold import TSNE
import umap

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# PCA
pca = PCA(n_components=2, random_state=42)
X_pca = pca.fit_transform(X_scaled)

# UMAP
reducer = umap.UMAP(n_components=2, random_state=42, n_neighbors=15, min_dist=0.1)
X_umap = reducer.fit_transform(X_scaled)

# t-SNE
tsne = TSNE(n_components=2, random_state=42, perplexity=30, max_iter=1000)
X_tsne = tsne.fit_transform(X_scaled)

各手法の特徴

手法	分散説明率	特徴
PCA	PC1: 44.3%, PC2: 19.0%（合計 63.2%）	線形射影。グローバル構造を保持。2 成分で分散の 63% を説明
UMAP	-	非線形。局所構造とグローバル構造の両方を保持。最も明確なクラスタ分離
t-SNE	-	非線形。局所構造に強い。クラスタの形状はパラメータ依存

3 手法すべてで悪性（赤）と良性（青）の明確な分離が確認できます。特に UMAP は 2 つのクラスタが最もコンパクトに分離しており、このデータセットの構造を最もよく捉えています。PCA の第 1 主成分だけで分散の 44.3% を説明しており、線形手法でも十分な判別力があることがわかります。

Step 4: 機械学習分類

4-1. モデルの訓練と評価

4 つの分類モデルを 80/20 の層化分割で訓練し、5-fold 層化交差検証で評価しました。

from sklearn.model_selection import train_test_split, StratifiedKFold, cross_val_score
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.svm import SVC
from xgboost import XGBClassifier

X_train, X_test, y_train, y_test = train_test_split(
    X_scaled, y, test_size=0.2, stratify=y, random_state=42
)

models = {
    "Logistic Regression": LogisticRegression(max_iter=5000, random_state=42),
    "Random Forest": RandomForestClassifier(n_estimators=200, max_depth=10, random_state=42),
    "SVM (RBF)": SVC(kernel="rbf", probability=True, random_state=42),
    "XGBoost": XGBClassifier(n_estimators=200, max_depth=5, learning_rate=0.1, random_state=42),
}

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

性能比較

モデル	CV AUC (mean +/- std)	Test AUC	Accuracy	F1 Score	Precision	Recall
Logistic Regression	0.9956 +/- 0.0050	0.9954	0.9825	0.9861	0.9861	0.9861
Random Forest	0.9894 +/- 0.0086	0.9931	0.9561	0.9655	0.9589	0.9722
SVM (RBF)	0.9956 +/- 0.0048	0.9950	0.9825	0.9861	0.9861	0.9861
XGBoost	0.9934 +/- 0.0049	0.9934	0.9474	0.9589	0.9459	0.9722

最優秀モデル: Logistic Regression（Test ROC-AUC = 0.9954）

なぜ Logistic Regression が最も高い性能か？
このデータセットは 30 特徴量に対して 569 サンプルと比較的低次元であり、PCA の結果からもわかるように線形分離可能に近い構造を持っています。このような場合、正則化された線形モデルは過学習を避けつつ高い汎化性能を発揮します。Random Forest や XGBoost のような複雑なモデルは、より高次元・非線形なデータで真価を発揮します。

4-2. ROC 曲線

4 モデルすべてが AUC > 0.99 を達成しており、左上隅に密集しています。ランダム分類器（対角線）との差は非常に大きく、データセットの判別可能性の高さを反映しています。

4-3. 混同行列

Logistic Regression と SVM は同じ性能（Accuracy = 0.983）を示し、テストセット 114 サンプル中 2 サンプルのみを誤分類しています。臨床応用の観点では、**偽陰性（悪性を良性と誤判定）**を最小化することが重要であり、Recall の高いモデルが望ましいです。

4-4. 特徴量重要度（Random Forest）

Random Forest の Gini 重要度から、worst concave points、worst perimeter、mean concave points が上位に位置しています。これは Step 2 の Mann-Whitney U 検定の結果とも一致しており、細胞核の凹みの深さと周囲長が悪性/良性の判別に最も重要な特徴量であることが、異なるアプローチから確認されました。

Step 5: 生存解析（合成データ）

5-1. 合成生存データの生成

実際の臨床データを模擬するため、200 患者の合成生存データを生成しました。

np.random.seed(42)
n_patients = 200
subtypes = np.random.choice(["ER+", "ER-", "HER2+"], size=n_patients, p=[0.5, 0.3, 0.2])

サブタイプ	患者数	イベント率	中央生存時間の設計
ER+（エストロゲン受容体陽性）	102	35%	scale=60 ヶ月（最良予後）
ER-（エストロゲン受容体陰性）	56	60%	scale=30 ヶ月（最悪予後）
HER2+（HER2 過発現）	42	45%	scale=45 ヶ月（中間予後）

合計イベント数: 83 / 200（41.5%）
打ち切り: 117 / 200（58.5%）
最大追跡期間: 120 ヶ月（10 年）

5-2. Kaplan-Meier 生存曲線

Log-rank 検定結果

比較	p-value	有意性
ER+ vs ER-	1.93e-08	***
ER+ vs HER2+	3.04e-01	ns
ER- vs HER2+	3.92e-04	***

ER-（赤）は ER+（青）および HER2+（緑）と比較して有意に予後不良であり、これは実際のがん臨床データの傾向と一致しています。ER+ と HER2+ の間には統計的有意差がありませんが、Kaplan-Meier 曲線を見ると ER+ のほうがやや良い傾向が見られます。

5-3. Cox 比例ハザードモデル

Cox PH モデルの結果

共変量	係数 (coef)	exp(coef) = HR	p-value	解釈
subtype_ER-	1.43	4.18	< 0.005	ER+ 基準で死亡リスク 4.18 倍
subtype_HER2+	0.31	1.36	0.29	統計的有意差なし
grade	0.11	1.11	0.48	統計的有意差なし
age	0.00	1.00	0.85	影響なし
tumor_size_cm	-0.02	0.98	0.67	影響なし

Concordance index = 0.70 — 中程度の予測能力
ER- サブタイプのみが統計的に有意な予後因子（HR = 4.18, p < 0.005）
age、tumor_size、grade は合成データのため有意な関連なし

合成データの限界: age、tumor_size、grade は生存時間と独立に生成しているため有意な係数が出ていません。実際の臨床データでは、これらの変数は生存時間と強く関連します。合成データはあくまで解析パイプラインの実証目的で使用しています。

実験から得られた知見

1. 線形モデルの強さ — 「シンプルが最強」の実例

4 モデルの比較で **Logistic Regression（AUC = 0.9954）**が最良という結果は、データサイエンスの重要な教訓を示しています。30 特徴量・569 サンプルという条件では、Random Forest や XGBoost のような複雑なモデルよりも、正則化された線形モデルが過学習を回避し、より高い汎化性能を発揮しました。モデルの複雑さとデータの複雑さを一致させることの重要性が確認できます。

2. 統計検定と機械学習の一致 — 解析の頑健性

Mann-Whitney U 検定で最も有意だった特徴量（worst concave points、worst perimeter）が、Random Forest の Gini 重要度でも上位に位置しました。異なるアプローチが同じ結論を支持することは、結果の信頼性を高める重要な要素です。

3. 次元削減の使い分け

手法	適したユースケース	本実験での知見
PCA	全体構造の把握、前処理	2 成分で分散の 63% を説明。線形分離可能性を示唆
UMAP	クラスタ発見、可視化	最も明確なクラスタ分離。大規模データにもスケーラブル
t-SNE	局所構造の可視化	perplexity パラメータに敏感。解釈に注意が必要

4. 生存解析パイプラインの汎用性

lifelines ライブラリを使った Kaplan-Meier + Cox PH の組み合わせは、合成データでも実データでも同じコードで実行可能です。ER- サブタイプの有意な予後不良（HR = 4.18）が正しく検出されたことで、パイプラインの妥当性が確認できました。

5. 再現性の確保

パイプライン全体で random_state=42 を統一し、matplotlib の Agg バックエンドを使用することで、環境に依存しない完全な再現性を実現しています。

まとめ

項目	内容
データ	scikit-learn Breast Cancer Wisconsin（569 samples x 30 features）
統計検定	Mann-Whitney U（全 30 特徴量が padj < 0.05）
次元削減	PCA（PC1+PC2 = 63.2%）、UMAP、t-SNE
最優秀モデル	Logistic Regression（Test AUC = 0.9954, Accuracy = 98.25%）
生存解析	合成データ 200 患者、ER- サブタイプ HR = 4.18 (p < 0.005)
図	10 publication-quality figures（300 DPI）
実行時間	約 30 秒

Claude Code を使うことで、がん遺伝子発現データの分類と生存解析の完全パイプラインを自然言語の指示だけで構築・実行できました。EDA からボルケーノプロット、4 モデルの比較、Kaplan-Meier 曲線、Cox 回帰まで、がんゲノミクス研究で必要とされる主要な解析を網羅しています。特に、統計検定と機械学習の結果が一致したことで、パイプラインの科学的妥当性が確認できた点は重要です。

参考資料

scikit-learn: Breast Cancer Wisconsin Dataset — データセットの公式ドキュメント
scikit-learn User Guide: Classification — 分類アルゴリズムのガイド
UMAP Documentation — Uniform Manifold Approximation and Projection
lifelines Documentation — Survival Analysis in Python
XGBoost Documentation — Gradient Boosting Framework
Street, W.N. et al. "Nuclear feature extraction for breast tumor diagnosis." IS&T/SPIE Symposium on Electronic Imaging 1993 — 原論文
McInnes, L. et al. "UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction." arXiv:1802.03426 (2018)
Davidson-Pilon, C. "lifelines: survival analysis in Python." JOSS 4(40), 1317 (2019)
Cox, D.R. "Regression Models and Life-Tables." Journal of the Royal Statistical Society B 34(2), 187-220 (1972)

ソースコード

本記事で使用した完全な解析スクリプトは以下の通りです。

cancer_classification.py（完全な解析パイプライン）

#!/usr/bin/env python3
"""
Experiment 03: Cancer Gene Expression Classification & Survival Analysis
========================================================================
Comprehensive cancer genomics analysis pipeline using scikit-learn's
breast cancer dataset and synthetic survival data.

Author: Claude Code (Exp-03)
Date: 2026-02-11
"""

import warnings
warnings.filterwarnings("ignore")

import matplotlib
matplotlib.use("Agg")

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import seaborn as sns
from pathlib import Path

# scikit-learn
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split, StratifiedKFold, cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier
from sklearn.svm import SVC
from sklearn.metrics import (
    roc_curve, auc, confusion_matrix, classification_report,
    accuracy_score, f1_score, precision_score, recall_score
)
from sklearn.decomposition import PCA

# Statistical tests
from scipy.stats import mannwhitneyu
from scipy.cluster.hierarchy import linkage

# Dimensionality reduction
import umap

# t-SNE
from sklearn.manifold import TSNE

# Survival analysis
from lifelines import KaplanMeierFitter, CoxPHFitter
from lifelines.statistics import logrank_test

# Volcano plot text adjustment
from adjustText import adjust_text

# XGBoost
try:
    from xgboost import XGBClassifier
    HAS_XGBOOST = True
except ImportError:
    HAS_XGBOOST = False

# ============================================================
# Configuration
# ============================================================
RANDOM_STATE = 42
np.random.seed(RANDOM_STATE)
DPI = 300
FIG_DIR = Path("/home/nahisaho/05_Experiment/ScientificSkills/Exp-03/figures")
FIG_DIR.mkdir(parents=True, exist_ok=True)

# Plot style
plt.rcParams.update({
    "figure.dpi": DPI,
    "savefig.dpi": DPI,
    "font.size": 10,
    "axes.titlesize": 12,
    "axes.labelsize": 10,
    "figure.facecolor": "white",
})

print("=" * 70)
print("Experiment 03: Cancer Gene Expression Classification & Survival Analysis")
print("=" * 70)

# ============================================================
# Step 1: Load Data & Exploratory Data Analysis
# ============================================================
print("\n[Step 1] Loading data and performing EDA...")

data = load_breast_cancer()
X = pd.DataFrame(data.data, columns=data.feature_names)
y = pd.Series(data.target, name="diagnosis")
target_names = data.target_names

print(f"  Dataset shape: {X.shape}")
print(f"  Classes: {dict(zip(target_names, np.bincount(y)))}")
print(f"  Features: {X.columns.tolist()[:5]}... ({len(X.columns)} total)")

# --- Fig 1: Feature distributions by diagnosis ---
top_features_for_dist = [
    "mean radius", "mean texture", "mean perimeter", "mean area",
    "mean smoothness", "mean compactness"
]

fig, axes = plt.subplots(2, 3, figsize=(14, 8))
fig.suptitle("Feature Distributions by Diagnosis", fontsize=14, fontweight="bold")

for idx, feat in enumerate(top_features_for_dist):
    ax = axes[idx // 3, idx % 3]
    for label_idx, label_name in enumerate(target_names):
        mask = y == label_idx
        ax.hist(
            X.loc[mask, feat], bins=30, alpha=0.6,
            label=label_name.capitalize(),
            color=["#e74c3c", "#3498db"][label_idx],
            edgecolor="white", linewidth=0.5
        )
    ax.set_title(feat, fontweight="bold")
    ax.set_xlabel("Value")
    ax.set_ylabel("Count")
    ax.legend(fontsize=8)

plt.tight_layout()
plt.savefig(FIG_DIR / "Fig1_feature_distributions.png", bbox_inches="tight")
plt.close()
print("  -> Fig1_feature_distributions.png saved")

# --- Fig 2: Correlation heatmap ---
top_var_features = X.var().nlargest(15).index.tolist()
corr_matrix = X[top_var_features].corr()

fig, ax = plt.subplots(figsize=(12, 10))
mask = np.triu(np.ones_like(corr_matrix, dtype=bool))
cmap = sns.diverging_palette(250, 10, as_cmap=True)
sns.heatmap(
    corr_matrix, mask=mask, cmap=cmap, center=0,
    annot=True, fmt=".2f", linewidths=0.5,
    square=True, ax=ax, annot_kws={"size": 7},
    cbar_kws={"shrink": 0.8, "label": "Pearson r"}
)
ax.set_title("Correlation Heatmap (Top 15 Features by Variance)",
             fontsize=13, fontweight="bold", pad=15)
plt.tight_layout()
plt.savefig(FIG_DIR / "Fig2_correlation_heatmap.png", bbox_inches="tight")
plt.close()
print("  -> Fig2_correlation_heatmap.png saved")

# ============================================================
# Step 2: Feature Selection & Differential Expression Analysis
# ============================================================
print("\n[Step 2] Differential expression analysis...")

malignant_mask = y == 0
benign_mask = y == 1

results = []
for feat in X.columns:
    mal_vals = X.loc[malignant_mask, feat]
    ben_vals = X.loc[benign_mask, feat]
    stat, pval = mannwhitneyu(mal_vals, ben_vals, alternative="two-sided")
    mean_mal = mal_vals.mean()
    mean_ben = ben_vals.mean()
    fc = mean_mal / max(mean_ben, 1e-10)
    log2fc = np.log2(fc) if fc > 0 else 0
    results.append({
        "feature": feat,
        "mean_malignant": mean_mal,
        "mean_benign": mean_ben,
        "log2_fold_change": log2fc,
        "p_value": pval,
        "neg_log10_pval": -np.log10(max(pval, 1e-300))
    })

de_df = pd.DataFrame(results).sort_values("p_value")
de_df["p_adjusted"] = np.minimum(de_df["p_value"] * len(de_df), 1.0)

print("\n  Top 10 discriminating features:")
print(de_df[["feature", "log2_fold_change", "p_value", "p_adjusted"]].head(10).to_string(index=False))

# --- Fig 3: Volcano Plot ---
fig, ax = plt.subplots(figsize=(10, 8))

de_df["significant"] = (de_df["p_adjusted"] < 0.05) & (de_df["log2_fold_change"].abs() > 0.5)
de_df["direction"] = "NS"
de_df.loc[(de_df["p_adjusted"] < 0.05) & (de_df["log2_fold_change"] > 0.5), "direction"] = "Up in Malignant"
de_df.loc[(de_df["p_adjusted"] < 0.05) & (de_df["log2_fold_change"] < -0.5), "direction"] = "Down in Malignant"

colors = {"NS": "#bdc3c7", "Up in Malignant": "#e74c3c", "Down in Malignant": "#3498db"}
for direction, color in colors.items():
    mask = de_df["direction"] == direction
    ax.scatter(
        de_df.loc[mask, "log2_fold_change"],
        de_df.loc[mask, "neg_log10_pval"],
        c=color, alpha=0.7, s=60, label=direction, edgecolors="white", linewidth=0.5
    )

texts = []
top_to_label = de_df.nlargest(10, "neg_log10_pval")
for _, row in top_to_label.iterrows():
    short_name = row["feature"].replace("mean ", "").replace("worst ", "w.")
    texts.append(
        ax.text(row["log2_fold_change"], row["neg_log10_pval"],
                short_name, fontsize=7, ha="center")
    )
adjust_text(texts, ax=ax, arrowprops=dict(arrowstyle="-", color="gray", lw=0.5))

ax.axhline(-np.log10(0.05 / len(de_df)), color="gray", ls="--", lw=0.8, label="Bonferroni threshold")
ax.axvline(-0.5, color="gray", ls=":", lw=0.8)
ax.axvline(0.5, color="gray", ls=":", lw=0.8)
ax.set_xlabel("log$_2$(Fold Change) [Malignant / Benign]", fontsize=11)
ax.set_ylabel("-log$_{10}$(p-value)", fontsize=11)
ax.set_title("Volcano Plot: Differential Feature Expression", fontsize=13, fontweight="bold")
ax.legend(loc="upper left", fontsize=9)
plt.tight_layout()
plt.savefig(FIG_DIR / "Fig3_volcano_plot.png", bbox_inches="tight")
plt.close()
print("  -> Fig3_volcano_plot.png saved")

# ============================================================
# Step 3: Dimensionality Reduction
# ============================================================
print("\n[Step 3] Dimensionality reduction (PCA, UMAP, t-SNE)...")

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

pca = PCA(n_components=2, random_state=RANDOM_STATE)
X_pca = pca.fit_transform(X_scaled)
pca_full = PCA(random_state=RANDOM_STATE)
pca_full.fit(X_scaled)

print(f"  PCA: PC1 variance = {pca.explained_variance_ratio_[0]:.3f}, "
      f"PC2 variance = {pca.explained_variance_ratio_[1]:.3f}")

reducer = umap.UMAP(n_components=2, random_state=RANDOM_STATE, n_neighbors=15, min_dist=0.1)
X_umap = reducer.fit_transform(X_scaled)

tsne = TSNE(n_components=2, random_state=RANDOM_STATE, perplexity=30, max_iter=1000)
X_tsne = tsne.fit_transform(X_scaled)

# --- Fig 4: PCA + UMAP + t-SNE ---
fig, axes = plt.subplots(1, 3, figsize=(18, 5.5))
color_map = {0: "#e74c3c", 1: "#3498db"}
colors_arr = [color_map[yi] for yi in y]

axes[0].scatter(X_pca[:, 0], X_pca[:, 1], c=colors_arr, alpha=0.6, s=20, edgecolors="white", linewidth=0.3)
axes[0].set_xlabel(f"PC1 ({pca.explained_variance_ratio_[0]*100:.1f}%)")
axes[0].set_ylabel(f"PC2 ({pca.explained_variance_ratio_[1]*100:.1f}%)")
axes[0].set_title("PCA", fontweight="bold")

axes[1].scatter(X_umap[:, 0], X_umap[:, 1], c=colors_arr, alpha=0.6, s=20, edgecolors="white", linewidth=0.3)
axes[1].set_xlabel("UMAP1")
axes[1].set_ylabel("UMAP2")
axes[1].set_title("UMAP", fontweight="bold")

scatter = axes[2].scatter(X_tsne[:, 0], X_tsne[:, 1], c=colors_arr, alpha=0.6, s=20, edgecolors="white", linewidth=0.3)
axes[2].set_xlabel("t-SNE1")
axes[2].set_ylabel("t-SNE2")
axes[2].set_title("t-SNE", fontweight="bold")

from matplotlib.patches import Patch
legend_elements = [
    Patch(facecolor="#e74c3c", label="Malignant"),
    Patch(facecolor="#3498db", label="Benign")
]
fig.legend(handles=legend_elements, loc="upper center", ncol=2, fontsize=11,
           bbox_to_anchor=(0.5, 1.02))

plt.suptitle("Dimensionality Reduction Comparison", fontsize=14, fontweight="bold", y=1.06)
plt.tight_layout()
plt.savefig(FIG_DIR / "Fig4_pca_umap.png", bbox_inches="tight")
plt.close()
print("  -> Fig4_pca_umap.png saved")

# ============================================================
# Step 4: Machine Learning Classification
# ============================================================
print("\n[Step 4] Training classification models...")

X_train, X_test, y_train, y_test = train_test_split(
    X_scaled, y, test_size=0.2, stratify=y, random_state=RANDOM_STATE
)

print(f"  Train: {X_train.shape[0]} samples, Test: {X_test.shape[0]} samples")

models = {
    "Logistic Regression": LogisticRegression(max_iter=5000, random_state=RANDOM_STATE, C=1.0),
    "Random Forest": RandomForestClassifier(n_estimators=200, random_state=RANDOM_STATE, max_depth=10),
    "SVM (RBF)": SVC(kernel="rbf", probability=True, random_state=RANDOM_STATE, C=1.0),
}

if HAS_XGBOOST:
    models["XGBoost"] = XGBClassifier(
        n_estimators=200, max_depth=5, learning_rate=0.1,
        random_state=RANDOM_STATE, use_label_encoder=False,
        eval_metric="logloss", verbosity=0
    )
else:
    models["GradientBoosting"] = GradientBoostingClassifier(
        n_estimators=200, max_depth=5, learning_rate=0.1, random_state=RANDOM_STATE
    )

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=RANDOM_STATE)
model_results = {}

for name, model in models.items():
    cv_scores = cross_val_score(model, X_train, y_train, cv=cv, scoring="roc_auc")
    model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    y_proba = model.predict_proba(X_test)[:, 1]

    acc = accuracy_score(y_test, y_pred)
    f1 = f1_score(y_test, y_pred)
    prec = precision_score(y_test, y_pred)
    rec = recall_score(y_test, y_pred)
    fpr, tpr, _ = roc_curve(y_test, y_proba)
    roc_auc_val = auc(fpr, tpr)

    model_results[name] = {
        "model": model, "cv_auc_mean": cv_scores.mean(), "cv_auc_std": cv_scores.std(),
        "test_accuracy": acc, "test_f1": f1, "test_precision": prec, "test_recall": rec,
        "test_roc_auc": roc_auc_val, "fpr": fpr, "tpr": tpr, "y_pred": y_pred, "y_proba": y_proba,
    }

    print(f"  {name:25s} | CV AUC: {cv_scores.mean():.4f} +/- {cv_scores.std():.4f} | "
          f"Test AUC: {roc_auc_val:.4f} | Acc: {acc:.4f}")

# --- Fig 5: ROC Curves ---
fig, ax = plt.subplots(figsize=(8, 7))
line_styles = ["-", "--", "-.", ":"]
colors_roc = ["#e74c3c", "#3498db", "#2ecc71", "#f39c12"]

for idx, (name, res) in enumerate(model_results.items()):
    ax.plot(res["fpr"], res["tpr"], color=colors_roc[idx], lw=2, ls=line_styles[idx],
            label=f'{name} (AUC = {res["test_roc_auc"]:.3f})')

ax.plot([0, 1], [0, 1], "k--", lw=1, alpha=0.5, label="Random (AUC = 0.500)")
ax.set_xlabel("False Positive Rate", fontsize=11)
ax.set_ylabel("True Positive Rate", fontsize=11)
ax.set_title("ROC Curves: Model Comparison", fontsize=13, fontweight="bold")
ax.legend(loc="lower right", fontsize=9)
ax.set_xlim([-0.02, 1.02]); ax.set_ylim([-0.02, 1.02])
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(FIG_DIR / "Fig5_roc_curves.png", bbox_inches="tight")
plt.close()
print("  -> Fig5_roc_curves.png saved")

# --- Fig 6: Confusion Matrices ---
n_models = len(model_results)
fig, axes = plt.subplots(1, n_models, figsize=(4 * n_models, 4))
if n_models == 1: axes = [axes]

for idx, (name, res) in enumerate(model_results.items()):
    cm = confusion_matrix(y_test, res["y_pred"])
    sns.heatmap(cm, annot=True, fmt="d", cmap="Blues",
                xticklabels=["Malignant", "Benign"], yticklabels=["Malignant", "Benign"],
                ax=axes[idx], cbar=False)
    axes[idx].set_title(f"{name}\n(Acc: {res['test_accuracy']:.3f})", fontsize=10, fontweight="bold")
    axes[idx].set_ylabel("True Label" if idx == 0 else "")
    axes[idx].set_xlabel("Predicted Label")

plt.suptitle("Confusion Matrices", fontsize=14, fontweight="bold", y=1.02)
plt.tight_layout()
plt.savefig(FIG_DIR / "Fig6_confusion_matrices.png", bbox_inches="tight")
plt.close()
print("  -> Fig6_confusion_matrices.png saved")

# --- Fig 7: Feature Importance (Random Forest) ---
rf_model = model_results["Random Forest"]["model"]
importances = rf_model.feature_importances_
feature_imp_df = pd.DataFrame({
    "feature": data.feature_names, "importance": importances
}).sort_values("importance", ascending=True)

fig, ax = plt.subplots(figsize=(10, 8))
top_n = 15
top_features = feature_imp_df.tail(top_n)
colors_imp = plt.cm.RdYlBu_r(np.linspace(0.2, 0.8, top_n))
ax.barh(range(top_n), top_features["importance"], color=colors_imp, edgecolor="white", linewidth=0.5)
ax.set_yticks(range(top_n))
ax.set_yticklabels(top_features["feature"], fontsize=9)
ax.set_xlabel("Feature Importance (Gini)", fontsize=11)
ax.set_title("Random Forest: Top 15 Feature Importances", fontsize=13, fontweight="bold")
ax.grid(axis="x", alpha=0.3)
plt.tight_layout()
plt.savefig(FIG_DIR / "Fig7_feature_importance.png", bbox_inches="tight")
plt.close()
print("  -> Fig7_feature_importance.png saved")

# ============================================================
# Step 5: Survival Analysis (Synthetic Data)
# ============================================================
print("\n[Step 5] Survival analysis with synthetic data...")

np.random.seed(RANDOM_STATE)
n_patients = 200
subtypes = np.random.choice(["ER+", "ER-", "HER2+"], size=n_patients, p=[0.5, 0.3, 0.2])

survival_data = []
for i in range(n_patients):
    subtype = subtypes[i]
    if subtype == "ER+":
        time = np.random.exponential(scale=60) + np.random.uniform(6, 12)
        event_prob = 0.35
    elif subtype == "ER-":
        time = np.random.exponential(scale=30) + np.random.uniform(3, 6)
        event_prob = 0.60
    else:
        time = np.random.exponential(scale=45) + np.random.uniform(4, 8)
        event_prob = 0.45

    event = np.random.binomial(1, event_prob)
    time = min(time, 120)
    time = max(time, 1)

    age = np.clip(np.random.normal(55, 12), 25, 85)
    tumor_size = np.clip(np.random.exponential(scale=2.5) + 0.5, 0.5, 15)
    grade = np.random.choice([1, 2, 3], p=[0.2, 0.5, 0.3])

    survival_data.append({
        "patient_id": i + 1, "subtype": subtype, "time_months": round(time, 1),
        "event": event, "age": round(age, 1), "tumor_size_cm": round(tumor_size, 1), "grade": grade,
    })

surv_df = pd.DataFrame(survival_data)
print(f"  Survival dataset: {surv_df.shape}")
print(f"  Events: {surv_df['event'].sum()} / {n_patients} ({surv_df['event'].mean()*100:.1f}%)")

# --- Fig 8: Kaplan-Meier Curves ---
fig, ax = plt.subplots(figsize=(10, 7))
kmf = KaplanMeierFitter()
subtype_colors = {"ER+": "#3498db", "ER-": "#e74c3c", "HER2+": "#2ecc71"}

for subtype, color in subtype_colors.items():
    mask = surv_df["subtype"] == subtype
    kmf.fit(surv_df.loc[mask, "time_months"], surv_df.loc[mask, "event"], label=subtype)
    kmf.plot_survival_function(ax=ax, color=color, ci_show=True, linewidth=2)

subtypes_list = list(subtype_colors.keys())
print("\n  Log-rank test results:")
for i in range(len(subtypes_list)):
    for j in range(i + 1, len(subtypes_list)):
        s1, s2 = subtypes_list[i], subtypes_list[j]
        m1, m2 = surv_df["subtype"] == s1, surv_df["subtype"] == s2
        lr_result = logrank_test(
            surv_df.loc[m1, "time_months"], surv_df.loc[m2, "time_months"],
            surv_df.loc[m1, "event"], surv_df.loc[m2, "event"]
        )
        sig = "***" if lr_result.p_value < 0.001 else ("**" if lr_result.p_value < 0.01 else ("*" if lr_result.p_value < 0.05 else "ns"))
        print(f"    {s1} vs {s2}: p = {lr_result.p_value:.4e} {sig}")

ax.set_xlabel("Time (months)", fontsize=11)
ax.set_ylabel("Survival Probability", fontsize=11)
ax.set_title("Kaplan-Meier Survival Curves by Subtype", fontsize=13, fontweight="bold")
ax.legend(fontsize=11, loc="lower left")
ax.set_xlim([0, 125]); ax.set_ylim([0, 1.05])
ax.grid(True, alpha=0.3)
ax.text(0.98, 0.02, f"n = {n_patients} patients\nEvents: {surv_df['event'].sum()}",
        transform=ax.transAxes, fontsize=9, va="bottom", ha="right",
        bbox=dict(boxstyle="round,pad=0.3", facecolor="wheat", alpha=0.5))
plt.tight_layout()
plt.savefig(FIG_DIR / "Fig8_kaplan_meier.png", bbox_inches="tight")
plt.close()
print("  -> Fig8_kaplan_meier.png saved")

# --- Cox PH Model ---
print("\n  Fitting Cox PH model...")
cox_df = surv_df.copy()
cox_df = pd.get_dummies(cox_df, columns=["subtype"], drop_first=True, dtype=float)
cox_df = cox_df.drop(columns=["patient_id"])

cph = CoxPHFitter()
cph.fit(cox_df, duration_col="time_months", event_col="event")
print(cph.print_summary())

# --- Fig 9: Cox Forest Plot ---
fig, ax = plt.subplots(figsize=(10, 6))
summary = cph.summary
coefs = summary["coef"]
ci_lower = summary["coef lower 95%"]
ci_upper = summary["coef upper 95%"]
p_values = summary["p"]

sorted_idx = coefs.abs().sort_values().index
coefs = coefs[sorted_idx]
ci_lower = ci_lower[sorted_idx]
ci_upper = ci_upper[sorted_idx]
p_values = p_values[sorted_idx]

y_pos = range(len(coefs))
colors_cox = ["#e74c3c" if p < 0.05 else "#95a5a6" for p in p_values]

ax.errorbar(coefs, y_pos, xerr=[coefs - ci_lower, ci_upper - coefs],
            fmt="o", color="black", ecolor="gray", elinewidth=1.5, capsize=3, markersize=6)
for i, (c, col) in enumerate(zip(coefs, colors_cox)):
    ax.scatter(c, i, color=col, s=80, zorder=5, edgecolors="black", linewidth=0.5)

ax.axvline(0, color="black", ls="--", lw=0.8)
ax.set_yticks(list(y_pos))
ax.set_yticklabels([name.replace("subtype_", "Subtype: ") for name in coefs.index], fontsize=9)
ax.set_xlabel("log(Hazard Ratio) [Cox PH Coefficient]", fontsize=11)
ax.set_title("Cox Proportional Hazards: Forest Plot", fontsize=13, fontweight="bold")
ax.grid(axis="x", alpha=0.3)

for i, (pv, cu) in enumerate(zip(p_values, ci_upper)):
    sig_txt = f"p={pv:.3f}" if pv >= 0.001 else f"p={pv:.1e}"
    ax.text(ci_upper.max() + 0.05, i, sig_txt, fontsize=8, va="center")

plt.tight_layout()
plt.savefig(FIG_DIR / "Fig9_cox_forest_plot.png", bbox_inches="tight")
plt.close()
print("  -> Fig9_cox_forest_plot.png saved")

# ============================================================
# Step 6: Results Summary
# ============================================================
print("\n[Step 6] Results summary...")

fig, axes = plt.subplots(1, 2, figsize=(14, 6))
metrics_names = ["CV AUC", "Test AUC", "Accuracy", "F1", "Precision", "Recall"]
model_names = list(model_results.keys())
short_names = [n.replace("Logistic Regression", "LogReg").replace("Random Forest", "RF")
               .replace("SVM (RBF)", "SVM").replace("GradientBoosting", "GBM") for n in model_names]

metrics_data = np.array([
    [res["cv_auc_mean"], res["test_roc_auc"], res["test_accuracy"],
     res["test_f1"], res["test_precision"], res["test_recall"]]
    for res in model_results.values()
])

x = np.arange(len(metrics_names))
width = 0.18
colors_bar = ["#e74c3c", "#3498db", "#2ecc71", "#f39c12"]

for i, (sname, color) in enumerate(zip(short_names, colors_bar)):
    axes[0].bar(x + i * width, metrics_data[i], width, label=sname, color=color, edgecolor="white")

axes[0].set_xticks(x + width * (len(model_names) - 1) / 2)
axes[0].set_xticklabels(metrics_names, fontsize=9)
axes[0].set_ylabel("Score", fontsize=11)
axes[0].set_title("Classification Metrics Comparison", fontsize=12, fontweight="bold")
axes[0].legend(fontsize=9)
axes[0].set_ylim([0.85, 1.02])
axes[0].grid(axis="y", alpha=0.3)

angles = np.linspace(0, 2 * np.pi, len(metrics_names), endpoint=False).tolist()
angles += angles[:1]

axes[1] = fig.add_subplot(122, polar=True)
for i, (sname, color) in enumerate(zip(short_names, colors_bar)):
    values = metrics_data[i].tolist()
    values += values[:1]
    axes[1].plot(angles, values, "o-", linewidth=2, color=color, label=sname, markersize=4)
    axes[1].fill(angles, values, alpha=0.1, color=color)

axes[1].set_xticks(angles[:-1])
axes[1].set_xticklabels(metrics_names, fontsize=8)
axes[1].set_ylim([0.85, 1.02])
axes[1].set_title("Radar: Model Performance", fontsize=12, fontweight="bold", pad=20)
axes[1].legend(loc="lower right", fontsize=8, bbox_to_anchor=(1.3, -0.05))

plt.tight_layout()
plt.savefig(FIG_DIR / "Fig10_model_comparison.png", bbox_inches="tight")
plt.close()
print("  -> Fig10_model_comparison.png saved")

# --- Performance table ---
print("\n" + "=" * 90)
print("PERFORMANCE COMPARISON TABLE")
print("=" * 90)
perf_table = pd.DataFrame({
    "Model": model_names,
    "CV AUC (mean +/- std)": [f'{res["cv_auc_mean"]:.4f} +/- {res["cv_auc_std"]:.4f}' for res in model_results.values()],
    "Test AUC": [f'{res["test_roc_auc"]:.4f}' for res in model_results.values()],
    "Accuracy": [f'{res["test_accuracy"]:.4f}' for res in model_results.values()],
    "F1 Score": [f'{res["test_f1"]:.4f}' for res in model_results.values()],
    "Precision": [f'{res["test_precision"]:.4f}' for res in model_results.values()],
    "Recall": [f'{res["test_recall"]:.4f}' for res in model_results.values()],
})
print(perf_table.to_string(index=False))

best_model_name = max(model_results, key=lambda k: model_results[k]["test_roc_auc"])
best_auc = model_results[best_model_name]["test_roc_auc"]
print(f"\nBest Model: {best_model_name} (Test ROC-AUC = {best_auc:.4f})")

# Save results
results_dir = Path("/home/nahisaho/05_Experiment/ScientificSkills/Exp-03/results")
results_dir.mkdir(parents=True, exist_ok=True)

perf_table.to_csv(results_dir / "model_performance.csv", index=False)
de_df.to_csv(results_dir / "differential_expression.csv", index=False)
surv_df.to_csv(results_dir / "survival_data.csv", index=False)

import json
summary_stats = {
    "dataset_samples": X.shape[0], "dataset_features": X.shape[1],
    "n_malignant": int((y == 0).sum()), "n_benign": int((y == 1).sum()),
    "pca_var_pc1": float(pca.explained_variance_ratio_[0]),
    "pca_var_pc2": float(pca.explained_variance_ratio_[1]),
    "best_model": best_model_name, "best_test_auc": float(best_auc),
    "survival_patients": n_patients, "survival_events": int(surv_df["event"].sum()),
}
with open(results_dir / "summary_stats.json", "w") as f:
    json.dump(summary_stats, f, indent=2)

print("\n" + "=" * 70)
print("Experiment 03 completed successfully!")
print("=" * 70)