{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.12.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceType":"competition","sourceId":29653,"databundleVersionId":2420395,"isSourceIdPinned":false}],"dockerImageVersionId":31328,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# 🧠 Radiogenomics Analytics Framework\n## MRI → Radiomics → Genomic Prediction (MGMT Methylation)\n\n> **Dataset:** RSNA-MICCAI Brain Tumor Radiogenomic Classification  \n> **Task:** Predict MGMT promoter methylation status from multi-modal MRI  \n> **Framework:** End-to-end radiogenomics pipeline for research & thesis use  \n\n---\n\n### 🔬 What is Radiogenomics?\nRadiogenomics is an interdisciplinary field that studies the **statistical associations between medical imaging features (radiomics)** and **genomic/molecular characteristics** of tumors. Instead of invasive biopsies, we extract quantitative imaging biomarkers and learn to predict molecular status non-invasively.\n\n**Clinical relevance:**  \n- MGMT (O6-methylguanine-DNA methyltransferase) promoter methylation is a key prognostic and predictive biomarker in glioblastoma (GBM)\n- MGMT-methylated tumors respond better to temozolomide chemotherapy\n- Standard detection requires invasive tissue biopsy — imaging-based prediction is highly valuable\n\n### 📋 Pipeline Overview\n```\nDICOM MRI (4 modalities)\n        ↓\n  Preprocessing (normalize, resize)\n        ↓\n  ROI Approximation (center crop + intensity threshold)\n        ↓\n  Radiomics Feature Extraction (intensity, texture, shape)\n        ↓\n  Feature Selection (correlation, ANOVA, LASSO, RF importance)\n        ↓\n  Radiogenomic Association Analysis\n        ↓\n  Multi-modality Fusion + Modeling\n        ↓\n  Evaluation & Interpretation\n```","metadata":{}},{"cell_type":"markdown","source":"---\n## A. 🛠️ Setup & Installation\n\nWe install required libraries and configure the environment. All libraries are available on Kaggle.","metadata":{}},{"cell_type":"code","source":"!ls /kaggle/input/competitions","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:04:00.783690Z","iopub.execute_input":"2026-03-29T08:04:00.784036Z","iopub.status.idle":"2026-03-29T08:04:00.909290Z","shell.execute_reply.started":"2026-03-29T08:04:00.784003Z","shell.execute_reply":"2026-03-29T08:04:00.908044Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION A: SETUP\n# ─────────────────────────────────────────────────────────────\n\n# Install pydicom if not available (usually pre-installed on Kaggle)\nimport subprocess\nsubprocess.run(['pip', 'install', 'pydicom', '-q'])\n\n# ── Core Libraries ──────────────────────────────────────────\nimport os\nimport warnings\nimport random\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport matplotlib.gridspec as gridspec\nimport seaborn as sns\nwarnings.filterwarnings('ignore')\n\n# ── Medical Imaging ──────────────────────────────────────────\nimport pydicom\nfrom pydicom.errors import InvalidDicomError\n\n# ── Image Processing ─────────────────────────────────────────\nfrom skimage import exposure, filters, measure\nfrom skimage.feature import graycomatrix, graycoprops\nfrom sklearn.preprocessing import StandardScaler\nfrom scipy import stats\n\n# ── Machine Learning ─────────────────────────────────────────\nfrom sklearn.linear_model import LogisticRegression, LassoCV\nfrom sklearn.ensemble import RandomForestClassifier\nfrom sklearn.feature_selection import SelectKBest, f_classif, mutual_info_classif\nfrom sklearn.model_selection import StratifiedKFold, cross_val_score, cross_validate\nfrom sklearn.metrics import (\n    roc_auc_score, accuracy_score,\n    confusion_matrix, ConfusionMatrixDisplay,\n    roc_curve, classification_report\n)\nfrom sklearn.pipeline import Pipeline\n\n# ── Reproducibility ──────────────────────────────────────────\nSEED = 42\nnp.random.seed(SEED)\nrandom.seed(SEED)\n\n# ── Styling ──────────────────────────────────────────────────\nplt.rcParams.update({\n    'figure.facecolor': '#0f0f1a',\n    'axes.facecolor': '#1a1a2e',\n    'axes.edgecolor': '#444466',\n    'axes.labelcolor': '#e0e0ff',\n    'xtick.color': '#b0b0cc',\n    'ytick.color': '#b0b0cc',\n    'text.color': '#e0e0ff',\n    'grid.color': '#2a2a4a',\n    'grid.linestyle': '--',\n    'grid.alpha': 0.5,\n    'font.family': 'DejaVu Sans',\n    'font.size': 11,\n})\nPALETTE = ['#7b68ee', '#00bcd4', '#ff6b6b', '#4caf50', '#ffa726']\n\n# ── Dataset Path ─────────────────────────────────────────────\nBASE_PATH = '/kaggle/input/competitions/rsna-miccai-brain-tumor-radiogenomic-classification'\nTRAIN_PATH = os.path.join(BASE_PATH, 'train')\nMODALITIES  = ['FLAIR', 'T1w', 'T1wCE', 'T2w']\n\nprint('✅ Setup complete!')\nprint(f'📂 Dataset path: {BASE_PATH}')\nprint(f'🧬 Modalities: {MODALITIES}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:05:55.684961Z","iopub.execute_input":"2026-03-29T08:05:55.685813Z","iopub.status.idle":"2026-03-29T08:05:59.308576Z","shell.execute_reply.started":"2026-03-29T08:05:55.685774Z","shell.execute_reply":"2026-03-29T08:05:59.307389Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## B. 📊 Data Understanding\n\nWe load the clinical labels, inspect class distribution, and visualize MRI scans across the four modalities.","metadata":{}},{"cell_type":"code","source":"print(BASE_PATH)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:06:00.120162Z","iopub.execute_input":"2026-03-29T08:06:00.120460Z","iopub.status.idle":"2026-03-29T08:06:00.126086Z","shell.execute_reply.started":"2026-03-29T08:06:00.120435Z","shell.execute_reply":"2026-03-29T08:06:00.125171Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION B.1: LOAD LABELS\n# ─────────────────────────────────────────────────────────────\n\ntrain_df = pd.read_csv(os.path.join(BASE_PATH, 'train_labels.csv'))\nprint('=== Training Labels ===' )\nprint(train_df.head(10))\nprint(f'\\nShape: {train_df.shape}')\nprint(f'\\nMGMT label counts:')\nprint(train_df['MGMT_value'].value_counts())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:06:01.976919Z","iopub.execute_input":"2026-03-29T08:06:01.977277Z","iopub.status.idle":"2026-03-29T08:06:02.023024Z","shell.execute_reply.started":"2026-03-29T08:06:01.977239Z","shell.execute_reply":"2026-03-29T08:06:02.022207Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION B.2: CLASS DISTRIBUTION\n# ─────────────────────────────────────────────────────────────\n\nfig, axes = plt.subplots(1, 2, figsize=(14, 5))\nfig.suptitle('MGMT Methylation Status — Class Distribution', fontsize=14, fontweight='bold', color='#e0e0ff')\n\n# Bar chart\ncounts = train_df['MGMT_value'].value_counts()\nbars = axes[0].bar(\n    ['MGMT Unmethylated (0)', 'MGMT Methylated (1)'],\n    counts.values,\n    color=PALETTE[:2], edgecolor='white', linewidth=0.5\n)\nfor bar, val in zip(bars, counts.values):\n    axes[0].text(bar.get_x() + bar.get_width()/2, bar.get_height() + 2,\n                 str(val), ha='center', va='bottom', fontweight='bold', color='white')\naxes[0].set_title('Absolute Count', color='#e0e0ff')\naxes[0].set_ylabel('Number of Patients')\naxes[0].grid(axis='y', alpha=0.3)\n\n# Pie chart\naxes[1].pie(\n    counts.values,\n    labels=['Unmethylated\\n(MGMT−)', 'Methylated\\n(MGMT+)'],\n    colors=PALETTE[:2],\n    autopct='%1.1f%%',\n    startangle=90,\n    textprops={'color': 'white', 'fontsize': 12}\n)\naxes[1].set_title('Proportion', color='#e0e0ff')\n\nplt.tight_layout()\nplt.savefig('class_distribution.png', dpi=150, bbox_inches='tight')\nplt.show()\n\nmethylated_pct = counts[1] / counts.sum() * 100\nprint(f'\\n📊 Class balance: {methylated_pct:.1f}% methylated')\nprint('💡 Note: Slight class imbalance — we will use stratified cross-validation')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:06:13.559620Z","iopub.execute_input":"2026-03-29T08:06:13.560401Z","iopub.status.idle":"2026-03-29T08:06:14.150436Z","shell.execute_reply.started":"2026-03-29T08:06:13.560366Z","shell.execute_reply":"2026-03-29T08:06:14.149517Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION B.3: MRI VISUALIZATION PER MODALITY\n# ─────────────────────────────────────────────────────────────\n\ndef load_dicom_slice(patient_id: str, modality: str, slice_idx: int = None) -> np.ndarray | None:\n    \"\"\"\n    Load a single DICOM slice from a given patient and modality.\n    Returns a normalized float32 array, or None if loading fails.\n    \"\"\"\n    folder = os.path.join(TRAIN_PATH, str(patient_id).zfill(5), modality)\n    if not os.path.exists(folder):\n        return None\n    dcm_files = sorted([f for f in os.listdir(folder) if f.endswith('.dcm')])\n    if not dcm_files:\n        return None\n    # Default: use middle slice for best anatomical content\n    idx = slice_idx if slice_idx is not None else len(dcm_files) // 2\n    idx = min(idx, len(dcm_files) - 1)\n    try:\n        dcm = pydicom.dcmread(os.path.join(folder, dcm_files[idx]))\n        img = dcm.pixel_array.astype(np.float32)\n        # Min-max normalization to [0, 1]\n        img_min, img_max = img.min(), img.max()\n        if img_max > img_min:\n            img = (img - img_min) / (img_max - img_min)\n        return img\n    except Exception:\n        return None\n\n\n# ── Visualize 2 patients × 4 modalities ──────────────────────\nsample_ids = train_df['BraTS21ID'].sample(2, random_state=SEED).values\nmodality_cmaps = {'FLAIR': 'magma', 'T1w': 'gray', 'T1wCE': 'hot', 'T2w': 'viridis'}\n\nfig, axes = plt.subplots(len(sample_ids), len(MODALITIES), figsize=(18, 5 * len(sample_ids)))\nfig.suptitle('Multi-Modal MRI Visualization\\n(Each row = 1 patient, each column = 1 MRI modality)',\n             fontsize=14, fontweight='bold', color='#e0e0ff', y=1.01)\n\nmodality_descriptions = {\n    'FLAIR': 'Fluid Attenuated\\nInversion Recovery',\n    'T1w':   'T1-weighted\\n(anatomy)',\n    'T1wCE': 'T1w + Contrast\\n(tumor enhancement)',\n    'T2w':   'T2-weighted\\n(edema / CSF)'\n}\n\nfor row_i, pid in enumerate(sample_ids):\n    label = train_df.loc[train_df['BraTS21ID'] == pid, 'MGMT_value'].values[0]\n    label_str = 'Methylated ✓' if label == 1 else 'Unmethylated ✗'\n    for col_j, mod in enumerate(MODALITIES):\n        ax = axes[row_i][col_j]\n        img = load_dicom_slice(pid, mod)\n        if img is not None:\n            ax.imshow(img, cmap=modality_cmaps[mod], aspect='equal')\n            ax.set_title(\n                f'{mod}\\n{modality_descriptions[mod]}',\n                color='#e0e0ff', fontsize=10, pad=4\n            )\n        else:\n            ax.text(0.5, 0.5, 'No DICOM', ha='center', va='center',\n                    transform=ax.transAxes, color='#888')\n        if col_j == 0:\n            ax.set_ylabel(f'Patient {pid}\\nMGMT: {label_str}',\n                         color='#ffa726' if label == 1 else '#7b68ee', fontsize=10)\n        ax.axis('off')\n\nplt.tight_layout()\nplt.savefig('mri_modalities.png', dpi=150, bbox_inches='tight')\nplt.show()\nprint('\\n📌 Key differences between modalities:')\nprint('  • FLAIR: Suppresses CSF; highlights perilesional edema & infiltration')\nprint('  • T1w:   Standard anatomy; low signal in edema')\nprint('  • T1wCE: Contrast enhances blood-brain-barrier disruption (active tumor core)')\nprint('  • T2w:   High signal in water/edema; peritumoral zone well visible')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:06:22.074702Z","iopub.execute_input":"2026-03-29T08:06:22.075066Z","iopub.status.idle":"2026-03-29T08:06:25.160188Z","shell.execute_reply.started":"2026-03-29T08:06:22.075035Z","shell.execute_reply":"2026-03-29T08:06:25.159134Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## C. 🔧 Image Preprocessing\n\nWe build a robust preprocessing pipeline: load DICOM → normalize intensity → resize → stack modalities.\n\n**Why preprocessing matters:**\n- DICOM pixel values vary across scanners/protocols — normalization ensures comparability\n- Consistent spatial dimensions are required for feature extraction\n- Multi-modality stacking enables cross-modal feature fusion","metadata":{}},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION C: IMAGE PREPROCESSING\n# ─────────────────────────────────────────────────────────────\n\nimport cv2  # Available on Kaggle\n\nTARGET_SIZE = (128, 128)  # Balance between speed and spatial detail\n\ndef preprocess_slice(img: np.ndarray, target_size: tuple = TARGET_SIZE) -> np.ndarray:\n    \"\"\"\n    Full preprocessing chain for a single MRI slice:\n    1. Resize to standard dimensions\n    2. Z-score normalization (zero mean, unit variance)\n    3. Clip to ±3σ to reduce outlier influence\n    4. Rescale to [0, 1]\n    \"\"\"\n    # 1. Resize\n    img_resized = cv2.resize(img, target_size, interpolation=cv2.INTER_AREA)\n    # 2. Z-score normalization\n    mean, std = img_resized.mean(), img_resized.std()\n    if std > 0:\n        img_z = (img_resized - mean) / std\n        # 3. Clip at ±3σ\n        img_z = np.clip(img_z, -3, 3)\n        # 4. Rescale to [0, 1]\n        img_norm = (img_z + 3) / 6.0\n    else:\n        img_norm = np.zeros_like(img_resized)\n    return img_norm.astype(np.float32)\n\n\ndef load_patient_volume(\n    patient_id: str,\n    modalities: list = MODALITIES,\n    n_slices: int = 5\n) -> dict:\n    \"\"\"\n    Load and preprocess multiple slices from all modalities for one patient.\n    \n    Args:\n        patient_id: BraTS21 patient ID\n        modalities: list of MRI modality names\n        n_slices: number of central slices to aggregate\n    \n    Returns:\n        dict mapping modality → averaged preprocessed slice (H×W)\n    \"\"\"\n    volume = {}\n    for mod in modalities:\n        folder = os.path.join(TRAIN_PATH, str(patient_id).zfill(5), mod)\n        if not os.path.exists(folder):\n            volume[mod] = None\n            continue\n        dcm_files = sorted([f for f in os.listdir(folder) if f.endswith('.dcm')])\n        if not dcm_files:\n            volume[mod] = None\n            continue\n        # Central slices contain the most diagnostically relevant tumor anatomy\n        mid = len(dcm_files) // 2\n        half = n_slices // 2\n        slice_indices = range(max(0, mid - half), min(len(dcm_files), mid + half + 1))\n        slices = []\n        for idx in slice_indices:\n            try:\n                dcm = pydicom.dcmread(os.path.join(folder, dcm_files[idx]))\n                raw = dcm.pixel_array.astype(np.float32)\n                slices.append(preprocess_slice(raw))\n            except Exception:\n                continue\n        volume[mod] = np.mean(slices, axis=0) if slices else None\n    return volume\n\n\n# ── Preprocessing Visualization ───────────────────────────────\ntest_pid = train_df['BraTS21ID'].iloc[0]\nraw_img = load_dicom_slice(test_pid, 'FLAIR')\nproc_img = preprocess_slice(raw_img) if raw_img is not None else None\n\nif raw_img is not None and proc_img is not None:\n    fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n    fig.suptitle(f'Preprocessing Pipeline — Patient {test_pid} (FLAIR)',\n                 fontsize=13, fontweight='bold', color='#e0e0ff')\n\n    axes[0].imshow(raw_img, cmap='magma')\n    axes[0].set_title('1. Raw DICOM\\n(min-max [0,1])', color='#e0e0ff')\n    axes[0].axis('off')\n\n    axes[1].imshow(proc_img, cmap='magma')\n    axes[1].set_title('2. Z-score + Clip + Resize\\n(128×128)', color='#e0e0ff')\n    axes[1].axis('off')\n\n    # Histogram comparison\n    axes[2].hist(raw_img.flatten(), bins=60, alpha=0.6, color=PALETTE[0],\n                 label='Raw', density=True)\n    axes[2].hist(proc_img.flatten(), bins=60, alpha=0.6, color=PALETTE[1],\n                 label='Preprocessed', density=True)\n    axes[2].set_title('3. Intensity Histogram', color='#e0e0ff')\n    axes[2].legend()\n    axes[2].set_xlabel('Pixel Intensity')\n    axes[2].set_ylabel('Density')\n    axes[2].grid(True, alpha=0.3)\n\n    plt.tight_layout()\n    plt.savefig('preprocessing.png', dpi=150, bbox_inches='tight')\n    plt.show()\n\nprint('✅ Preprocessing pipeline ready')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:06:32.242020Z","iopub.execute_input":"2026-03-29T08:06:32.242340Z","iopub.status.idle":"2026-03-29T08:06:33.716042Z","shell.execute_reply.started":"2026-03-29T08:06:32.242316Z","shell.execute_reply":"2026-03-29T08:06:33.715104Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## D. 🎯 Tumor Region Approximation (Pseudo-ROI)\n\n**Real radiomics requires expert tumor segmentation** (e.g., BraTS annotations) to compute features only within the lesion volume.\n\n**Our pragmatic approach** uses two complementary strategies:\n1. **Center crop** — 60% central region where tumors predominantly occur\n2. **Intensity thresholding** — isolate bright hyperintense regions (FLAIR) that likely represent tumor/edema\n\n⚠️ **Limitations:**\n- May include normal brain tissue, reducing feature specificity\n- Cannot differentiate tumor core from edema without proper segmentation\n- Real pipelines use tools like FSL, FreeSurfer, or deep learning segmenters (nnU-Net)","metadata":{}},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION D: TUMOR REGION APPROXIMATION\n# ─────────────────────────────────────────────────────────────\n\ndef get_center_crop(img: np.ndarray, crop_frac: float = 0.6) -> np.ndarray:\n    \"\"\"\n    Extract central crop of the image.\n    Brain tumors in GBM patients predominantly occur in the central/peri-ventricular region.\n    \"\"\"\n    h, w = img.shape\n    ch, cw = int(h * crop_frac), int(w * crop_frac)\n    top  = (h - ch) // 2\n    left = (w - cw) // 2\n    return img[top:top+ch, left:left+cw]\n\n\ndef get_intensity_mask(img: np.ndarray, percentile: float = 75.0) -> np.ndarray:\n    \"\"\"\n    Binary mask: pixels brighter than the given percentile.\n    In FLAIR, tumor/edema appear hyperintense (bright) relative to normal brain.\n    Returns mask of same spatial dimensions as input.\n    \"\"\"\n    threshold = np.percentile(img, percentile)\n    return (img > threshold).astype(np.float32)\n\n\ndef apply_roi(\n    img: np.ndarray,\n    modality: str = 'FLAIR',\n    use_intensity: bool = True\n) -> np.ndarray:\n    \"\"\"\n    Combined ROI: center crop followed by intensity mask (modality-aware).\n    \n    - FLAIR, T2w: apply intensity threshold (hyperintense tumor/edema)\n    - T1wCE:      apply intensity threshold (enhancing tumor)\n    - T1w:        center crop only (less hyperintense contrast)\n    \"\"\"\n    cropped = get_center_crop(img)\n    if use_intensity and modality in ['FLAIR', 'T2w', 'T1wCE']:\n        percentile = 70 if modality == 'T1wCE' else 75\n        mask = get_intensity_mask(cropped, percentile)\n        return cropped * mask\n    return cropped\n\n\n# ── Visualize ROI approximation ───────────────────────────────\nif proc_img is not None:\n    center_crop = get_center_crop(proc_img)\n    intensity_mask = get_intensity_mask(proc_img)\n    roi_img = apply_roi(proc_img, modality='FLAIR')\n\n    fig, axes = plt.subplots(1, 4, figsize=(18, 4))\n    fig.suptitle('Tumor Region Approximation Strategies (FLAIR)',\n                 fontsize=13, fontweight='bold', color='#e0e0ff')\n\n    titles = ['Original\\n(preprocessed)', 'Center Crop\\n(60%)',\n              'Intensity Mask\\n(>75th pct)', 'Combined ROI\\n(crop × mask)']\n    imgs = [proc_img, center_crop, proc_img * intensity_mask, roi_img]\n\n    for ax, title, im in zip(axes, titles, imgs):\n        ax.imshow(im, cmap='magma')\n        ax.set_title(title, color='#e0e0ff', fontsize=10)\n        ax.axis('off')\n\n    plt.tight_layout()\n    plt.savefig('roi_approximation.png', dpi=150, bbox_inches='tight')\n    plt.show()\n\nprint('⚠️  ROI Limitation reminder:')\nprint('   This is a spatial heuristic. For publication-quality radiomics,')\nprint('   use provided segmentation masks or automated segmentation (nnU-Net).')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:06:41.735804Z","iopub.execute_input":"2026-03-29T08:06:41.736160Z","iopub.status.idle":"2026-03-29T08:06:42.335535Z","shell.execute_reply.started":"2026-03-29T08:06:41.736130Z","shell.execute_reply":"2026-03-29T08:06:42.334574Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## E. 🔬 Radiomics Feature Extraction\n\nWe extract three categories of quantitative imaging features from each patient:\n\n| Category | Features | Biological Relevance |\n|---|---|---|\n| **Intensity** | Mean, SD, skewness, kurtosis, energy, entropy | Tumor cellularity, necrosis, heterogeneity |\n| **Texture (GLCM)** | Contrast, correlation, homogeneity, dissimilarity | Tissue heterogeneity, architectural complexity |\n| **Shape/Morphology** | Area, perimeter, solidity, eccentricity | Tumor invasiveness, geometry |\n\nWe extract these for **each of the 4 modalities**, creating a rich multi-modal feature matrix.","metadata":{}},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION E: RADIOMICS FEATURE EXTRACTION\n# ─────────────────────────────────────────────────────────────\n\ndef extract_intensity_features(roi: np.ndarray, prefix: str = '') -> dict:\n    \"\"\"\n    First-order intensity statistics.\n    \n    - Mean:     average brightness → tissue density proxy\n    - Std:      intensity spread → intratumoral heterogeneity\n    - Skewness: distribution asymmetry → necrosis/enhancement patterns\n    - Kurtosis: tail behavior → outlier intensity regions\n    - Energy:   sum of squared intensities → overall signal strength\n    - Entropy:  information content → textural complexity\n    - P25/P75:  quartile intensities → robust range measure\n    - IQR:      interquartile range → spread without outliers\n    \"\"\"\n    flat = roi.flatten()\n    flat = flat[flat > 0]  # Exclude background (masked regions)\n    if len(flat) < 10:\n        return {f'{prefix}_{k}': np.nan for k in [\n            'mean','std','skewness','kurtosis','energy','entropy','p25','p75','iqr']}\n    p25, p75 = np.percentile(flat, [25, 75])\n    hist, _ = np.histogram(flat, bins=64, density=True)\n    hist = hist[hist > 0]\n    return {\n        f'{prefix}_mean':     float(np.mean(flat)),\n        f'{prefix}_std':      float(np.std(flat)),\n        f'{prefix}_skewness': float(stats.skew(flat)),\n        f'{prefix}_kurtosis': float(stats.kurtosis(flat)),\n        f'{prefix}_energy':   float(np.sum(flat ** 2)),\n        f'{prefix}_entropy':  float(-np.sum(hist * np.log2(hist + 1e-10))),\n        f'{prefix}_p25':      float(p25),\n        f'{prefix}_p75':      float(p75),\n        f'{prefix}_iqr':      float(p75 - p25),\n    }\n\n\ndef extract_glcm_features(roi: np.ndarray, prefix: str = '') -> dict:\n    \"\"\"\n    Gray-Level Co-occurrence Matrix (GLCM) texture features.\n    GLCM captures spatial relationships between pixel intensity pairs.\n    \n    - Contrast:      local intensity variation → edge sharpness\n    - Correlation:   linear dependency → structural regularity\n    - Homogeneity:   closeness to diagonal → texture smoothness\n    - Dissimilarity: weighted contrast → local texture variation\n    - Energy (ASM):  sum of squared elements → textural uniformity\n    \"\"\"\n    # Convert to 8-bit for GLCM (requires integer-valued input)\n    img_uint8 = (roi * 255).astype(np.uint8)\n    if img_uint8.max() == 0:\n        return {f'{prefix}_glcm_{k}': np.nan for k in [\n            'contrast','correlation','homogeneity','dissimilarity','energy']}\n    # Compute GLCM at 4 orientations → rotation-invariant features via mean\n    glcm = graycomatrix(\n        img_uint8, distances=[1, 3],\n        angles=[0, np.pi/4, np.pi/2, 3*np.pi/4],\n        levels=256, symmetric=True, normed=True\n    )\n    props = ['contrast', 'correlation', 'homogeneity', 'dissimilarity', 'energy']\n    return {\n        f'{prefix}_glcm_{prop}': float(graycoprops(glcm, prop).mean())\n        for prop in props\n    }\n\n\ndef extract_shape_features(roi: np.ndarray, prefix: str = '') -> dict:\n    \"\"\"\n    Morphological / shape proxy features from the ROI mask.\n    These capture the geometric properties of the hyperintense region,\n    which proxies the gross tumor shape.\n    \n    - Area fraction:  proportion of bright pixels → tumor burden proxy\n    - Solidity:       convex hull fill → convexity/invasion proxy\n    - Eccentricity:   elongation → shape irregularity\n    - Perimeter:      boundary length → surface complexity\n    - Compactness:    4π×area/perimeter² → circularity\n    \"\"\"\n    binary = (roi > roi.mean()).astype(np.uint8)\n    total_pixels = binary.size\n    area_frac = float(binary.sum() / total_pixels)\n    regions = measure.regionprops(binary)\n    if not regions:\n        return {\n            f'{prefix}_area_frac':   area_frac,\n            f'{prefix}_solidity':    np.nan,\n            f'{prefix}_eccentricity':np.nan,\n            f'{prefix}_perimeter':   np.nan,\n            f'{prefix}_compactness': np.nan,\n        }\n    r = max(regions, key=lambda x: x.area)  # Largest connected region\n    perimeter = r.perimeter if r.perimeter > 0 else 1\n    compactness = (4 * np.pi * r.area) / (perimeter ** 2)\n    return {\n        f'{prefix}_area_frac':    area_frac,\n        f'{prefix}_solidity':     float(r.solidity),\n        f'{prefix}_eccentricity': float(r.eccentricity),\n        f'{prefix}_perimeter':    float(perimeter),\n        f'{prefix}_compactness':  float(compactness),\n    }\n\n\ndef extract_all_features(patient_id: str) -> dict:\n    \"\"\"\n    Master feature extraction function for one patient.\n    Loads all modalities, applies ROI approximation, extracts all feature categories.\n    Returns a flat dictionary of all features.\n    \"\"\"\n    volume = load_patient_volume(patient_id)\n    all_features = {'BraTS21ID': patient_id}\n    for mod in MODALITIES:\n        img = volume.get(mod)\n        if img is None:\n            # Fill with NaN for missing modalities\n            dummy_prefix = mod\n            all_features.update({f'{dummy_prefix}_mean': np.nan})\n            continue\n        roi = apply_roi(img, modality=mod)\n        prefix = mod\n        all_features.update(extract_intensity_features(roi, prefix))\n        all_features.update(extract_glcm_features(roi, prefix))\n        all_features.update(extract_shape_features(roi, prefix))\n    return all_features\n\n\nprint('🔬 Feature extraction functions defined.')\nprint(f'   Features per modality: ~{9 + 5 + 5} = 19')\nprint(f'   Total features (4 modalities): ~{4 * 19} = 76')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:06:51.098498Z","iopub.execute_input":"2026-03-29T08:06:51.099187Z","iopub.status.idle":"2026-03-29T08:06:51.122609Z","shell.execute_reply.started":"2026-03-29T08:06:51.099148Z","shell.execute_reply":"2026-03-29T08:06:51.121302Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# FEATURE EXTRACTION: Run on N=50 patient subset\n# ─────────────────────────────────────────────────────────────\n# Using 50 patients keeps runtime manageable (<10 min on Kaggle)\n# For a full study, process all ~585 training patients\n\nN_PATIENTS = 50  # Increase to None for full dataset\n\n# Stratified sample to maintain class balance\nsubset_df = train_df.groupby('MGMT_value', group_keys=False).apply(\n    lambda x: x.sample(min(N_PATIENTS // 2, len(x)), random_state=SEED)\n).reset_index(drop=True)\n\nprint(f'📊 Using {len(subset_df)} patients ({subset_df[\"MGMT_value\"].value_counts().to_dict()})')\nprint('⏳ Extracting radiomics features... (may take 5-10 minutes)')\n\nfeature_records = []\nfailed = []\nfor i, row in subset_df.iterrows():\n    pid = row['BraTS21ID']\n    try:\n        feats = extract_all_features(pid)\n        feats['MGMT_value'] = row['MGMT_value']\n        feature_records.append(feats)\n        if (len(feature_records)) % 10 == 0:\n            print(f'  → Processed {len(feature_records)}/{len(subset_df)} patients')\n    except Exception as e:\n        failed.append(pid)\n\nfeat_df = pd.DataFrame(feature_records)\nprint(f'\\n✅ Feature matrix shape: {feat_df.shape}')\nprint(f'   Failed: {len(failed)} patients')\nfeat_df.head(3)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:06:57.605444Z","iopub.execute_input":"2026-03-29T08:06:57.605810Z","iopub.status.idle":"2026-03-29T08:07:22.015917Z","shell.execute_reply.started":"2026-03-29T08:06:57.605778Z","shell.execute_reply":"2026-03-29T08:07:22.015047Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Feature matrix overview ───────────────────────────────────\n\n# Separate feature columns from metadata\nmeta_cols = ['BraTS21ID', 'MGMT_value']\nfeature_cols = [c for c in feat_df.columns if c not in meta_cols]\n\nX_raw = feat_df[feature_cols].copy()\ny     = feat_df['MGMT_value'].values\n\nprint(f'Feature matrix: {X_raw.shape[0]} patients × {X_raw.shape[1]} features')\nprint(f'Missing values per feature (top 10):')\nprint(X_raw.isna().sum().sort_values(ascending=False).head(10))\n\n# Impute NaN with median (robust to outliers)\nX_raw = X_raw.fillna(X_raw.median())\nprint(f'\\n✅ After median imputation — NaN count: {X_raw.isna().sum().sum()}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:07:28.636581Z","iopub.execute_input":"2026-03-29T08:07:28.637149Z","iopub.status.idle":"2026-03-29T08:07:28.672419Z","shell.execute_reply.started":"2026-03-29T08:07:28.637114Z","shell.execute_reply":"2026-03-29T08:07:28.671407Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## F. 🔍 Feature Selection\n\nHigh-dimensional feature spaces cause overfitting. We apply a multi-stage selection pipeline:\n\n1. **Correlation filtering** — remove redundant features (r > 0.90)\n2. **ANOVA F-test** — univariate statistical test per feature vs MGMT\n3. **Mutual Information** — non-linear dependency measure\n4. **LASSO** — L1-regularized logistic regression (sparse selection)\n5. **Random Forest importance** — ensemble-based ranking","metadata":{}},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION F: FEATURE SELECTION\n# ─────────────────────────────────────────────────────────────\n\nfrom sklearn.preprocessing import StandardScaler\n\n# ── Step 1: Correlation Filtering ────────────────────────────\nprint('━━ Step 1: Correlation Filtering ━━')\ncorr_matrix = X_raw.corr().abs()\nupper_tri = corr_matrix.where(\n    np.triu(np.ones(corr_matrix.shape), k=1).astype(bool)\n)\n# Remove one feature from each highly correlated pair\nto_drop = [col for col in upper_tri.columns if any(upper_tri[col] > 0.90)]\nX_decorr = X_raw.drop(columns=to_drop)\nprint(f'  Original features: {X_raw.shape[1]}')\nprint(f'  Dropped (corr > 0.90): {len(to_drop)}')\nprint(f'  Remaining: {X_decorr.shape[1]}')\n\n# ── Step 2: Standardize ───────────────────────────────────────\nscaler = StandardScaler()\nX_scaled = scaler.fit_transform(X_decorr)\nfeature_names_decorr = X_decorr.columns.tolist()\n\n# ── Step 3: ANOVA F-test ──────────────────────────────────────\nprint('\\n━━ Step 2: ANOVA F-test ━━')\nanova_selector = SelectKBest(score_func=f_classif, k='all')\nanova_selector.fit(X_scaled, y)\nf_scores  = anova_selector.scores_\nf_pvalues = anova_selector.pvalues_\n\nanova_df = pd.DataFrame({\n    'feature': feature_names_decorr,\n    'f_score': f_scores,\n    'p_value': f_pvalues\n}).sort_values('f_score', ascending=False)\n\nsig_anova = anova_df[anova_df['p_value'] < 0.05]\nprint(f'  Significant features (p < 0.05): {len(sig_anova)}')\nprint(anova_df.head(10).to_string(index=False))\n\n# ── Step 4: Mutual Information ────────────────────────────────\nprint('\\n━━ Step 3: Mutual Information ━━')\nmi_scores = mutual_info_classif(X_scaled, y, random_state=SEED)\nmi_df = pd.DataFrame({\n    'feature': feature_names_decorr,\n    'mi_score': mi_scores\n}).sort_values('mi_score', ascending=False)\nprint(f'  Top MI features:')\nprint(mi_df.head(10).to_string(index=False))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:07:41.430569Z","iopub.execute_input":"2026-03-29T08:07:41.431671Z","iopub.status.idle":"2026-03-29T08:07:41.593381Z","shell.execute_reply.started":"2026-03-29T08:07:41.431605Z","shell.execute_reply":"2026-03-29T08:07:41.592305Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Step 5: LASSO Feature Selection ──────────────────────────\nprint('━━ Step 4: LASSO Regularization ━━')\nfrom sklearn.linear_model import LassoCV\n\nlasso = LassoCV(cv=5, random_state=SEED, max_iter=5000)\nlasso.fit(X_scaled, y)\nlasso_coefs = pd.Series(np.abs(lasso.coef_), index=feature_names_decorr)\nlasso_selected = lasso_coefs[lasso_coefs > 0].sort_values(ascending=False)\nprint(f'  Features selected by LASSO (coef ≠ 0): {len(lasso_selected)}')\nprint(lasso_selected.head(15))\n\n# ── Step 6: Random Forest Importance ─────────────────────────\nprint('\\n━━ Step 5: Random Forest Feature Importance ━━')\nrf_fs = RandomForestClassifier(n_estimators=200, random_state=SEED, n_jobs=-1)\nrf_fs.fit(X_scaled, y)\nrf_importances = pd.Series(rf_fs.feature_importances_, index=feature_names_decorr)\nrf_top = rf_importances.sort_values(ascending=False).head(20)\nprint(f'  Top 20 Random Forest features:')\nprint(rf_top)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:07:54.473406Z","iopub.execute_input":"2026-03-29T08:07:54.473782Z","iopub.status.idle":"2026-03-29T08:07:57.719300Z","shell.execute_reply.started":"2026-03-29T08:07:54.473752Z","shell.execute_reply":"2026-03-29T08:07:57.718157Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Consensus Feature Set ─────────────────────────────────────\n# Features that appear in top-20 by at least 2 of 3 methods: ANOVA, MI, RF\n\nTOP_K = 20\nanova_top = set(anova_df.head(TOP_K)['feature'].tolist())\nmi_top    = set(mi_df.head(TOP_K)['feature'].tolist())\nrf_top_set = set(rf_top.index.tolist())\nlasso_set  = set(lasso_selected.index.tolist())\n\n# Consensus: appear in ≥2 of 4 methods\nfrom collections import Counter\nvote_counter = Counter()\nfor feat_set in [anova_top, mi_top, rf_top_set, lasso_set]:\n    for f in feat_set:\n        vote_counter[f] += 1\n\nconsensus_features = [f for f, v in vote_counter.items() if v >= 2]\nprint(f'\\n📌 Consensus features (≥2 methods agree): {len(consensus_features)}')\nprint(consensus_features)\n\n# Final feature matrices\nX_selected = X_decorr[consensus_features].values\nX_selected = scaler.fit_transform(X_selected)  # Re-scale on final set\n\n# ── Feature Selection Summary Plot ───────────────────────────\nfig, axes = plt.subplots(1, 3, figsize=(20, 6))\nfig.suptitle('Feature Selection — Multi-Method Ranking',\n             fontsize=14, fontweight='bold', color='#e0e0ff')\n\n# ANOVA\ntop10 = anova_df.head(15)\naxes[0].barh(top10['feature'][::-1], top10['f_score'][::-1], color=PALETTE[0])\naxes[0].set_title('ANOVA F-score (Top 15)', color='#e0e0ff')\naxes[0].set_xlabel('F-score')\naxes[0].grid(axis='x', alpha=0.3)\n\n# Mutual Information\nmi_top15 = mi_df.head(15)\naxes[1].barh(mi_top15['feature'][::-1], mi_top15['mi_score'][::-1], color=PALETTE[1])\naxes[1].set_title('Mutual Information (Top 15)', color='#e0e0ff')\naxes[1].set_xlabel('MI Score')\naxes[1].grid(axis='x', alpha=0.3)\n\n# RF Importance\nrf_top15 = rf_importances.sort_values(ascending=False).head(15)\naxes[2].barh(rf_top15.index[::-1], rf_top15.values[::-1], color=PALETTE[2])\naxes[2].set_title('Random Forest Importance (Top 15)', color='#e0e0ff')\naxes[2].set_xlabel('Importance')\naxes[2].grid(axis='x', alpha=0.3)\n\nfor ax in axes:\n    ax.tick_params(labelsize=8)\n\nplt.tight_layout()\nplt.savefig('feature_selection.png', dpi=150, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:08:03.285690Z","iopub.execute_input":"2026-03-29T08:08:03.286721Z","iopub.status.idle":"2026-03-29T08:08:04.587877Z","shell.execute_reply.started":"2026-03-29T08:08:03.286685Z","shell.execute_reply":"2026-03-29T08:08:04.586909Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## G. 🧬 Radiogenomic Analysis (CRITICAL SECTION)\n\nThis section is the heart of radiogenomics: **linking imaging features to genomic/molecular status**.\n\nWe perform:\n1. **Statistical association testing** — which features correlate with MGMT status?\n2. **Multi-modality fusion analysis** — does combining MRI sequences improve prediction?\n3. **Feature interpretation** — which features are biologically plausible MGMT markers?","metadata":{}},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION G.1: STATISTICAL ASSOCIATION ANALYSIS\n# ─────────────────────────────────────────────────────────────\n\nprint('═══ RADIOGENOMIC ASSOCIATION ANALYSIS ═══')\nprint('Testing: which imaging features are statistically associated with MGMT methylation?\\n')\n\nresults = []\nfor feat in feature_names_decorr:\n    vals = X_decorr[feat].values\n    group0 = vals[y == 0]\n    group1 = vals[y == 1]\n    group0 = group0[~np.isnan(group0)]\n    group1 = group1[~np.isnan(group1)]\n    if len(group0) < 3 or len(group1) < 3:\n        continue\n    # Mann-Whitney U test (non-parametric, robust to non-normality)\n    stat, pval = stats.mannwhitneyu(group0, group1, alternative='two-sided')\n    # Effect size: rank-biserial correlation\n    n1, n2 = len(group0), len(group1)\n    effect_size = 1 - (2 * stat) / (n1 * n2)\n    results.append({\n        'feature':     feat,\n        'u_statistic': stat,\n        'p_value':     pval,\n        'effect_size': abs(effect_size),\n        'mean_unmeth': group0.mean(),\n        'mean_meth':   group1.mean(),\n        'delta_mean':  group1.mean() - group0.mean(),\n    })\n\nassoc_df = pd.DataFrame(results).sort_values('p_value')\n\n# FDR correction (Benjamini-Hochberg)\nfrom statsmodels.stats.multitest import multipletests\n_, fdr_pvals, _, _ = multipletests(assoc_df['p_value'].values, method='fdr_bh')\nassoc_df['fdr_pval'] = fdr_pvals\nassoc_df['significant'] = assoc_df['fdr_pval'] < 0.05\n\nprint('Top 15 features associated with MGMT methylation:')\ndisplay_cols = ['feature', 'p_value', 'fdr_pval', 'effect_size', 'mean_unmeth', 'mean_meth', 'significant']\nprint(assoc_df[display_cols].head(15).to_string(index=False, float_format='{:.4f}'.format))\n\nn_sig = assoc_df['significant'].sum()\nprint(f'\\n🔬 Significant after FDR correction (q < 0.05): {n_sig} features')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:08:14.380186Z","iopub.execute_input":"2026-03-29T08:08:14.381005Z","iopub.status.idle":"2026-03-29T08:08:14.501475Z","shell.execute_reply.started":"2026-03-29T08:08:14.380972Z","shell.execute_reply":"2026-03-29T08:08:14.500280Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Volcano Plot: Radiogenomic Association ────────────────────\n\nfig, axes = plt.subplots(1, 2, figsize=(18, 6))\nfig.suptitle('Radiogenomic Association: Imaging Features vs MGMT Methylation',\n             fontsize=14, fontweight='bold', color='#e0e0ff')\n\n# Volcano Plot\nax = axes[0]\nlog_p  = -np.log10(assoc_df['p_value'].clip(1e-10))\neffect = assoc_df['effect_size']\ncolors_pt = ['#ff6b6b' if sig else '#444466' for sig in assoc_df['significant']]\n\nax.scatter(effect, log_p, c=colors_pt, alpha=0.8, s=40, edgecolors='none')\nax.axhline(-np.log10(0.05), color='#ffa726', linestyle='--', alpha=0.7, label='p=0.05')\nax.set_xlabel('Effect Size (rank-biserial |r|)', color='#e0e0ff')\nax.set_ylabel('-log10(p-value)', color='#e0e0ff')\nax.set_title('Volcano Plot — Feature–MGMT Associations', color='#e0e0ff')\nax.legend()\nax.grid(True, alpha=0.3)\n\n# Annotate top 5 features\nfor _, row in assoc_df.head(5).iterrows():\n    feat_name = row['feature'].replace('_', '\\n')\n    xi = row['effect_size']\n    yi = -np.log10(max(row['p_value'], 1e-10))\n    ax.annotate(feat_name, (xi, yi), fontsize=7, color='#ffd700',\n                xytext=(xi+0.02, yi+0.1), arrowprops=dict(arrowstyle='->', color='#888'))\n\n# Box plots: top 5 significant features\nax2 = axes[1]\ntop5_feats = assoc_df.head(5)['feature'].tolist()\nplot_data_meth  = [X_decorr.loc[y == 1, f].dropna().values for f in top5_feats]\nplot_data_unmeth = [X_decorr.loc[y == 0, f].dropna().values for f in top5_feats]\n\npositions_0 = np.arange(len(top5_feats)) * 2\npositions_1 = positions_0 + 0.7\n\nbp0 = ax2.boxplot(plot_data_unmeth, positions=positions_0, widths=0.6,\n                  patch_artist=True,\n                  boxprops=dict(facecolor='#7b68ee', alpha=0.7),\n                  medianprops=dict(color='white', linewidth=2),\n                  whiskerprops=dict(color='#7b68ee'),\n                  capprops=dict(color='#7b68ee'),\n                  flierprops=dict(marker='o', color='#7b68ee', alpha=0.3, markersize=3))\nbp1 = ax2.boxplot(plot_data_meth, positions=positions_1, widths=0.6,\n                  patch_artist=True,\n                  boxprops=dict(facecolor='#ff6b6b', alpha=0.7),\n                  medianprops=dict(color='white', linewidth=2),\n                  whiskerprops=dict(color='#ff6b6b'),\n                  capprops=dict(color='#ff6b6b'),\n                  flierprops=dict(marker='o', color='#ff6b6b', alpha=0.3, markersize=3))\n\nax2.set_xticks(positions_0 + 0.35)\nax2.set_xticklabels([f.replace('_', '\\n') for f in top5_feats], fontsize=8)\nax2.set_title('Top 5 Features: MGMT− vs MGMT+', color='#e0e0ff')\nax2.set_ylabel('Feature Value (standardized)')\nax2.legend([bp0['boxes'][0], bp1['boxes'][0]], ['MGMT− (unmethylated)', 'MGMT+ (methylated)'],\n           loc='upper right')\nax2.grid(True, alpha=0.3)\n\nplt.tight_layout()\nplt.savefig('radiogenomic_association.png', dpi=150, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:08:20.349684Z","iopub.execute_input":"2026-03-29T08:08:20.350680Z","iopub.status.idle":"2026-03-29T08:08:21.825957Z","shell.execute_reply.started":"2026-03-29T08:08:20.350604Z","shell.execute_reply":"2026-03-29T08:08:21.825039Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION G.2: MULTI-MODALITY FUSION ANALYSIS\n# ─────────────────────────────────────────────────────────────\n# Compare predictive power: each modality alone vs all combined\n\nprint('═══ MULTI-MODALITY FUSION ANALYSIS ═══')\nprint('Hypothesis: combining FLAIR+T1w+T1wCE+T2w outperforms any single modality\\n')\n\ncv = StratifiedKFold(n_splits=5, shuffle=True, random_state=SEED)\n\ndef evaluate_modality_subset(modality_list: list, feat_df: pd.DataFrame, y: np.ndarray,\n                               cv=cv, label: str = '') -> dict:\n    \"\"\"\n    Evaluate a Logistic Regression model using features from specified modalities only.\n    Returns mean ROC-AUC, accuracy from 5-fold CV.\n    \"\"\"\n    mod_feats = [c for c in feat_df.columns\n                 if any(c.startswith(mod + '_') for mod in modality_list)]\n    if not mod_feats:\n        return {'label': label, 'auc': np.nan, 'acc': np.nan}\n    X_mod = feat_df[mod_feats].fillna(feat_df[mod_feats].median()).values\n    pipe = Pipeline([\n        ('scaler', StandardScaler()),\n        ('clf', LogisticRegression(max_iter=1000, random_state=SEED, C=0.1))\n    ])\n    scores = cross_validate(pipe, X_mod, y, cv=cv,\n                            scoring=['roc_auc', 'accuracy'], n_jobs=-1)\n    return {\n        'label':    label,\n        'auc':      scores['test_roc_auc'].mean(),\n        'auc_std':  scores['test_roc_auc'].std(),\n        'acc':      scores['test_accuracy'].mean(),\n        'acc_std':  scores['test_accuracy'].std(),\n        'n_feats':  len(mod_feats)\n    }\n\nfusion_results = []\nfor mod in MODALITIES:\n    res = evaluate_modality_subset([mod], X_decorr, y, label=f'{mod} only')\n    fusion_results.append(res)\n    print(f'  {mod:8s}: AUC = {res[\"auc\"]:.3f} ± {res[\"auc_std\"]:.3f}  |  ACC = {res[\"acc\"]:.3f}')\n\n# All modalities combined\nres_all = evaluate_modality_subset(MODALITIES, X_decorr, y, label='All Modalities (fusion)')\nfusion_results.append(res_all)\nprint(f'  {\"All fused\":8s}: AUC = {res_all[\"auc\"]:.3f} ± {res_all[\"auc_std\"]:.3f}  |  ACC = {res_all[\"acc\"]:.3f}')\n\nfusion_df = pd.DataFrame(fusion_results)\nprint(f'\\n🏆 Best single modality: {fusion_df.iloc[:-1].loc[fusion_df.iloc[:-1][\"auc\"].idxmax(), \"label\"]}')\nprint(f'🔀 Multi-modality fusion AUC: {res_all[\"auc\"]:.3f}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:08:28.558428Z","iopub.execute_input":"2026-03-29T08:08:28.558790Z","iopub.status.idle":"2026-03-29T08:08:31.231286Z","shell.execute_reply.started":"2026-03-29T08:08:28.558759Z","shell.execute_reply":"2026-03-29T08:08:31.230357Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Fusion Analysis Bar Chart ─────────────────────────────────\n\nfig, ax = plt.subplots(figsize=(12, 5))\nx_pos = np.arange(len(fusion_df))\nbar_colors = PALETTE[:4] + ['#ffd700']\nbars = ax.bar(x_pos, fusion_df['auc'], color=bar_colors, alpha=0.85,\n              yerr=fusion_df['auc_std'], capsize=5, edgecolor='white', linewidth=0.5)\nax.set_xticks(x_pos)\nax.set_xticklabels(fusion_df['label'], fontsize=11)\nax.set_ylabel('ROC-AUC (5-fold CV mean ± std)')\nax.set_title('Multi-Modality Fusion Analysis\\nSingle Modality vs All-Modality Fusion',\n             color='#e0e0ff', fontweight='bold')\nax.set_ylim(0.3, 1.0)\nax.axhline(0.5, color='#888', linestyle='--', alpha=0.6, label='Random baseline')\nfor bar, val, std in zip(bars, fusion_df['auc'], fusion_df['auc_std']):\n    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + std + 0.01,\n            f'{val:.3f}', ha='center', fontsize=10, fontweight='bold', color='white')\nax.legend()\nax.grid(axis='y', alpha=0.3)\nplt.tight_layout()\nplt.savefig('modality_fusion.png', dpi=150, bbox_inches='tight')\nplt.show()\n\nprint('\\n📌 Clinical interpretation:')\nprint('   FLAIR captures perilesional edema — closely linked to tumor infiltration')\nprint('   T1wCE reveals blood-brain barrier disruption — MGMT links to vascular abnormalities')\nprint('   Multi-modal fusion leverages complementary biological information')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:08:45.354230Z","iopub.execute_input":"2026-03-29T08:08:45.354564Z","iopub.status.idle":"2026-03-29T08:08:45.830074Z","shell.execute_reply.started":"2026-03-29T08:08:45.354535Z","shell.execute_reply":"2026-03-29T08:08:45.829039Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION G.3: FEATURE INTERPRETATION\n# Biological meaning of top radiogenomic features\n# ─────────────────────────────────────────────────────────────\n\nBIOLOGICAL_ANNOTATIONS = {\n    'mean':        'Average signal intensity → tissue cellularity, edema volume',\n    'entropy':     'Information heterogeneity → intratumoral heterogeneity (ITH), key MGMT correlate',\n    'skewness':    'Signal asymmetry → necrotic vs viable tumor fraction',\n    'kurtosis':    'Tail behavior → outlier enhancement zones, micro-necrosis',\n    'std':         'Signal variance → heterogeneous microenvironment',\n    'glcm_contrast':'Textural contrast → architectural irregularity, cellular density variation',\n    'glcm_correlation': 'Structural regularity → tissue organization (lower in aggressive GBM)',\n    'glcm_homogeneity': 'Textural smoothness → uniform tissue regions',\n    'glcm_energy': 'Textural uniformity → homogeneous tumor sub-regions',\n    'glcm_dissimilarity': 'Local heterogeneity → boundary irregularity',\n    'area_frac':   'Proportion of hyperintense region → tumor/edema burden',\n    'solidity':    'Convexity → tumor margin regularity (irregular in MGMT-unmethylated)',\n    'eccentricity':'Shape elongation → infiltrative growth pattern',\n    'compactness': 'Shape circularity → degree of tumor boundary regularity',\n    'energy':      'Total signal energy → absolute tumor signal load',\n    'iqr':         'Interquartile intensity range → robust heterogeneity measure',\n}\n\nprint('═══ BIOLOGICAL INTERPRETATION OF TOP FEATURES ═══\\n')\nfor feat in assoc_df.head(10)['feature']:\n    # Extract the core feature type from the prefixed name\n    for key, bio in BIOLOGICAL_ANNOTATIONS.items():\n        if feat.endswith(key):\n            mod = feat.split('_')[0]\n            pval = assoc_df.loc[assoc_df['feature'] == feat, 'p_value'].values[0]\n            eff  = assoc_df.loc[assoc_df['feature'] == feat, 'effect_size'].values[0]\n            print(f'📍 {feat}')\n            print(f'   Modality: {mod}  |  p={pval:.4f}  |  effect={eff:.3f}')\n            print(f'   Biology:  {bio}\\n')\n            break","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:08:54.250930Z","iopub.execute_input":"2026-03-29T08:08:54.251257Z","iopub.status.idle":"2026-03-29T08:08:54.267409Z","shell.execute_reply.started":"2026-03-29T08:08:54.251225Z","shell.execute_reply":"2026-03-29T08:08:54.266304Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## H. 🤖 Predictive Modeling\n\nWe compare:\n- **Logistic Regression** — interpretable linear model, good baseline\n- **Random Forest** — ensemble, captures non-linear interactions\n- **Before vs after feature selection** — quantify selection benefit\n\nAll evaluated with **5-fold stratified cross-validation** (no data leakage).","metadata":{}},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION H: MODELING\n# ─────────────────────────────────────────────────────────────\n\nfrom sklearn.dummy import DummyClassifier\n\ncv = StratifiedKFold(n_splits=5, shuffle=True, random_state=SEED)\n\n# ── Define models ─────────────────────────────────────────────\nmodels = {\n    'Dummy (baseline)': DummyClassifier(strategy='stratified', random_state=SEED),\n    'Logistic Regression': LogisticRegression(max_iter=2000, C=0.1, random_state=SEED),\n    'Random Forest':  RandomForestClassifier(n_estimators=300, max_depth=5,\n                                             random_state=SEED, n_jobs=-1),\n}\n\n# ── Feature configurations ────────────────────────────────────\nfeature_configs = {\n    'All Features (before selection)': X_scaled,  # Full decorrelated set\n    'Selected Features (after selection)': X_selected,  # Consensus set\n}\n\nresults_table = []\n\nprint('═══ MODEL EVALUATION ═══')\nprint(f'Cross-validation: {cv.n_splits}-fold stratified\\n')\n\nfor feat_name, X_feat in feature_configs.items():\n    print(f'\\n── Features: {feat_name} ({X_feat.shape[1]} features) ──')\n    for model_name, clf in models.items():\n        pipe = Pipeline([('clf', clf)])\n        cv_results = cross_validate(\n            pipe, X_feat, y, cv=cv,\n            scoring=['roc_auc', 'accuracy', 'f1'],\n            n_jobs=-1\n        )\n        auc  = cv_results['test_roc_auc'].mean()\n        auc_std = cv_results['test_roc_auc'].std()\n        acc  = cv_results['test_accuracy'].mean()\n        f1   = cv_results['test_f1'].mean()\n        results_table.append({\n            'Features': feat_name.split(' (')[0],\n            'Model': model_name,\n            'AUC': auc,\n            'AUC_std': auc_std,\n            'Accuracy': acc,\n            'F1': f1\n        })\n        print(f'  {model_name:30s}: AUC={auc:.3f}±{auc_std:.3f}  Acc={acc:.3f}  F1={f1:.3f}')\n\nresults_df = pd.DataFrame(results_table)\nprint('\\n✅ Modeling complete')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:09:03.557531Z","iopub.execute_input":"2026-03-29T08:09:03.557902Z","iopub.status.idle":"2026-03-29T08:09:07.263588Z","shell.execute_reply.started":"2026-03-29T08:09:03.557860Z","shell.execute_reply":"2026-03-29T08:09:07.262602Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Performance Comparison Plot ───────────────────────────────\n\nfig, axes = plt.subplots(1, 2, figsize=(16, 6))\nfig.suptitle('Model Comparison: Before vs After Feature Selection',\n             fontsize=14, fontweight='bold', color='#e0e0ff')\n\nfor ax_i, metric in enumerate(['AUC', 'Accuracy']):\n    ax = axes[ax_i]\n    pivot = results_df.pivot(index='Model', columns='Features', values=metric)\n    x = np.arange(len(pivot.index))\n    w = 0.35\n    cols = pivot.columns.tolist()\n    for j, col in enumerate(cols):\n        offset = (j - 0.5) * w\n        bars = ax.bar(x + offset, pivot[col], w, label=col,\n                     color=PALETTE[j], alpha=0.85, edgecolor='white', linewidth=0.5)\n        if metric == 'AUC' and col in results_df['Features'].values:\n            stds = results_df[(results_df['Features'] == col)]['AUC_std'].values\n            ax.errorbar(x + offset, pivot[col], yerr=stds, fmt='none',\n                        ecolor='white', capsize=4, alpha=0.7)\n    ax.set_xticks(x)\n    ax.set_xticklabels(pivot.index, rotation=15, ha='right', fontsize=9)\n    ax.set_ylabel(metric)\n    ax.set_title(metric, color='#e0e0ff')\n    ax.set_ylim(0.3, 1.05)\n    ax.axhline(0.5, color='#888', linestyle='--', alpha=0.4)\n    ax.legend(fontsize=8)\n    ax.grid(axis='y', alpha=0.3)\n\nplt.tight_layout()\nplt.savefig('model_comparison.png', dpi=150, bbox_inches='tight')\nplt.show()\n\nprint('\\n📋 Summary Table:')\nprint(results_df.to_string(index=False, float_format='{:.3f}'.format))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:09:15.319898Z","iopub.execute_input":"2026-03-29T08:09:15.320238Z","iopub.status.idle":"2026-03-29T08:09:16.116056Z","shell.execute_reply.started":"2026-03-29T08:09:15.320210Z","shell.execute_reply":"2026-03-29T08:09:16.114897Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## I. 📈 Evaluation\n\nDetailed evaluation of the best model: ROC curve, confusion matrix, and classification report.","metadata":{}},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION I: EVALUATION\n# ─────────────────────────────────────────────────────────────\n\n# Re-train best model on full data for visualization\n# (In practice, evaluation should stay within CV folds)\nbest_model = RandomForestClassifier(n_estimators=300, max_depth=5, random_state=SEED, n_jobs=-1)\nbest_model.fit(X_selected, y)\n\n# ── ROC Curve via CV ─────────────────────────────────────────\nfig, axes = plt.subplots(1, 3, figsize=(20, 6))\nfig.suptitle('Model Evaluation — Random Forest (Selected Features)',\n             fontsize=14, fontweight='bold', color='#e0e0ff')\n\n# ROC curves per fold\nax_roc = axes[0]\nmean_fpr = np.linspace(0, 1, 100)\ntprs, aucs = [], []\n\nfor fold, (train_idx, val_idx) in enumerate(cv.split(X_selected, y)):\n    clf_fold = RandomForestClassifier(n_estimators=300, max_depth=5, random_state=SEED, n_jobs=-1)\n    clf_fold.fit(X_selected[train_idx], y[train_idx])\n    y_prob = clf_fold.predict_proba(X_selected[val_idx])[:, 1]\n    fpr, tpr, _ = roc_curve(y[val_idx], y_prob)\n    auc_val = roc_auc_score(y[val_idx], y_prob)\n    aucs.append(auc_val)\n    interp_tpr = np.interp(mean_fpr, fpr, tpr)\n    interp_tpr[0] = 0.0\n    tprs.append(interp_tpr)\n    ax_roc.plot(fpr, tpr, alpha=0.3, color=PALETTE[fold % len(PALETTE)],\n                label=f'Fold {fold+1} (AUC={auc_val:.3f})')\n\nmean_tpr = np.mean(tprs, axis=0)\nmean_tpr[-1] = 1.0\nmean_auc = np.mean(aucs)\nstd_auc  = np.std(aucs)\nax_roc.plot(mean_fpr, mean_tpr, color='white', linewidth=2.5,\n            label=f'Mean ROC (AUC={mean_auc:.3f}±{std_auc:.3f})')\nax_roc.fill_between(mean_fpr,\n                    np.mean(tprs, axis=0) - np.std(tprs, axis=0),\n                    np.mean(tprs, axis=0) + np.std(tprs, axis=0),\n                    alpha=0.15, color='white')\nax_roc.plot([0,1],[0,1], 'r--', alpha=0.5, label='Random')\nax_roc.set_xlabel('False Positive Rate')\nax_roc.set_ylabel('True Positive Rate')\nax_roc.set_title('ROC Curves (5-Fold CV)', color='#e0e0ff')\nax_roc.legend(fontsize=7)\nax_roc.grid(True, alpha=0.3)\n\n# Confusion matrix (on full data for display only)\ny_pred = best_model.predict(X_selected)\ncm = confusion_matrix(y, y_pred)\ncm_disp = ConfusionMatrixDisplay(cm, display_labels=['MGMT−\\n(Unmethylated)', 'MGMT+\\n(Methylated)'])\ncm_disp.plot(ax=axes[1], colorbar=False, cmap='Blues')\naxes[1].set_title('Confusion Matrix\\n(full data, for display)', color='#e0e0ff')\naxes[1].set_facecolor('#1a1a2e')\n\n# AUC per fold bar chart\naxes[2].bar(range(1, 6), aucs, color=PALETTE[:5], alpha=0.85, edgecolor='white')\naxes[2].axhline(mean_auc, color='white', linewidth=2, linestyle='--', label=f'Mean={mean_auc:.3f}')\naxes[2].axhline(0.5, color='#ff6b6b', linewidth=1, linestyle=':', alpha=0.7, label='Random=0.5')\naxes[2].set_xlabel('CV Fold')\naxes[2].set_ylabel('ROC-AUC')\naxes[2].set_title('AUC per CV Fold', color='#e0e0ff')\naxes[2].set_ylim(0, 1)\naxes[2].legend()\naxes[2].grid(axis='y', alpha=0.3)\n\nplt.tight_layout()\nplt.savefig('evaluation.png', dpi=150, bbox_inches='tight')\nplt.show()\n\n# Classification report\nprint('\\n📋 Classification Report (full data):')\nprint(classification_report(y, y_pred, target_names=['MGMT− (0)', 'MGMT+ (1)']))\nprint(f'Mean CV AUC: {mean_auc:.3f} ± {std_auc:.3f}')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:09:19.988239Z","iopub.execute_input":"2026-03-29T08:09:19.988573Z","iopub.status.idle":"2026-03-29T08:09:25.316419Z","shell.execute_reply.started":"2026-03-29T08:09:19.988546Z","shell.execute_reply":"2026-03-29T08:09:25.315496Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"---\n## J. 🔍 Interpretation & Biological Relevance\n\nFeature importance analysis with biological annotation — connecting imaging patterns to molecular biology.","metadata":{}},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION J: INTERPRETATION\n# ─────────────────────────────────────────────────────────────\n\n# ── Feature Importance (Random Forest) ───────────────────────\nrf_importances_sel = pd.Series(\n    best_model.feature_importances_,\n    index=consensus_features\n).sort_values(ascending=True)\n\nfig = plt.figure(figsize=(18, 12))\ngs = gridspec.GridSpec(2, 2, figure=fig, hspace=0.4, wspace=0.35)\nfig.suptitle('Radiogenomics Feature Interpretation Dashboard',\n             fontsize=15, fontweight='bold', color='#e0e0ff', y=1.01)\n\n# ── Panel 1: Feature Importance ───────────────────────────────\nax1 = fig.add_subplot(gs[0, 0])\ncolors_bar = []\nfor f in rf_importances_sel.index:\n    if f.startswith('FLAIR'): colors_bar.append(PALETTE[0])\n    elif f.startswith('T1w'):  colors_bar.append(PALETTE[1])\n    elif f.startswith('T1wCE'):colors_bar.append(PALETTE[2])\n    else:                     colors_bar.append(PALETTE[3])\n\nbars = ax1.barh(range(len(rf_importances_sel)), rf_importances_sel.values,\n                color=colors_bar, alpha=0.85, edgecolor='none')\nax1.set_yticks(range(len(rf_importances_sel)))\nax1.set_yticklabels([f.replace('_glcm_', '\\nGLCM_').replace('_', '\\n')\n                     for f in rf_importances_sel.index], fontsize=7)\nax1.set_xlabel('Feature Importance')\nax1.set_title('Random Forest Feature Importance\\n(selected features only)', color='#e0e0ff', fontsize=10)\n# Legend\nlegend_patches = [\n    plt.Rectangle((0,0),1,1, fc=PALETTE[i], label=m)\n    for i, m in enumerate(MODALITIES)\n]\nax1.legend(handles=legend_patches, loc='lower right', fontsize=7)\nax1.grid(axis='x', alpha=0.3)\n\n# ── Panel 2: Correlation heatmap ──────────────────────────────\nax2 = fig.add_subplot(gs[0, 1])\nfeat_subset = consensus_features[:min(12, len(consensus_features))]\ncorr_subset = X_decorr[feat_subset].corr()\nmask = np.triu(np.ones_like(corr_subset, dtype=bool), k=1)\nsns.heatmap(\n    corr_subset, ax=ax2, mask=mask,\n    cmap='coolwarm', center=0, vmin=-1, vmax=1,\n    annot=True, fmt='.2f', annot_kws={'size': 6},\n    xticklabels=[f.replace('_', '\\n') for f in feat_subset],\n    yticklabels=[f.replace('_', '\\n') for f in feat_subset],\n    cbar_kws={'shrink': 0.7}\n)\nax2.set_title('Feature Correlation Matrix\\n(selected features)', color='#e0e0ff', fontsize=10)\nax2.tick_params(labelsize=6)\n\n# ── Panel 3: MGMT group comparison radar ──────────────────────\nax3 = fig.add_subplot(gs[1, 0])\ntop5 = rf_importances_sel.tail(5).index.tolist()\nfor label, color, ls in [(0, PALETTE[0], '-'), (1, PALETTE[2], '--')]:\n    group_means = X_decorr.loc[y == label, top5].mean()\n    ax3.plot(range(len(top5)), group_means.values, color=color, linewidth=2,\n             linestyle=ls, marker='o', markersize=6,\n             label=f'MGMT {\"−\" if label==0 else \"+\"} (n={sum(y==label)})')\n    ax3.fill_between(range(len(top5)), group_means.values, alpha=0.1, color=color)\nax3.set_xticks(range(len(top5)))\nax3.set_xticklabels([f.replace('_', '\\n') for f in top5], fontsize=8)\nax3.set_title('Top 5 Features: MGMT− vs MGMT+\\n(mean feature values)', color='#e0e0ff', fontsize=10)\nax3.set_ylabel('Mean Feature Value')\nax3.legend()\nax3.grid(True, alpha=0.3)\n\n# ── Panel 4: Modality contribution pie ───────────────────────\nax4 = fig.add_subplot(gs[1, 1])\nmod_importance = {}\nfor mod in MODALITIES:\n    mod_feats_sel = [f for f in consensus_features if f.startswith(mod + '_')]\n    mod_importance[mod] = sum(best_model.feature_importances_[consensus_features.index(f)]\n                               for f in mod_feats_sel if f in consensus_features)\n\nmods_sorted = sorted(mod_importance, key=mod_importance.get, reverse=True)\nax4.pie(\n    [mod_importance[m] for m in mods_sorted],\n    labels=mods_sorted,\n    colors=PALETTE[:4],\n    autopct='%1.1f%%',\n    startangle=90,\n    textprops={'color': 'white', 'fontsize': 11}\n)\nax4.set_title('Modality Contribution to\\nRandom Forest Prediction', color='#e0e0ff', fontsize=10)\n\nplt.savefig('interpretation_dashboard.png', dpi=150, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:09:56.742495Z","iopub.execute_input":"2026-03-29T08:09:56.743012Z","iopub.status.idle":"2026-03-29T08:09:59.761091Z","shell.execute_reply.started":"2026-03-29T08:09:56.742977Z","shell.execute_reply":"2026-03-29T08:09:59.760163Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# SECTION J.2: BIOLOGICAL RELEVANCE SUMMARY\n# ─────────────────────────────────────────────────────────────\n\nprint('═══════════════════════════════════════════════════════════')\nprint('  RADIOGENOMICS BIOLOGICAL INTERPRETATION SUMMARY')\nprint('═══════════════════════════════════════════════════════════\\n')\n\nprint('''\n🧬 MGMT METHYLATION — BIOLOGICAL BACKGROUND\n─────────────────────────────────────────────\nMGMT (O6-methylguanine-DNA methyltransferase) is a DNA repair enzyme.\nWhen its promoter is methylated (silenced), the repair mechanism is\ninactivated, making the tumor more sensitive to alkylating agents\n(e.g., temozolomide chemotherapy).\n\nMGMT methylation is present in ~40-50% of GBM patients and correlates\nwith improved overall survival under standard-of-care chemotherapy.\n\n🔬 WHY IMAGING CAN PREDICT MGMT STATUS\n────────────────────────────────────────\n1. ENTROPY / HETEROGENEITY:\n   MGMT-methylated GBMs tend to exhibit greater FLAIR signal heterogeneity,\n   reflecting complex mixtures of viable tumor, necrosis, and edema.\n   Literature shows entropy features are among the most discriminative\n   radiomics biomarkers for MGMT prediction.\n\n2. TEXTURE (GLCM CONTRAST, DISSIMILARITY):\n   MGMT-unmethylated tumors often show more homogeneous enhancement\n   (T1wCE) due to uniform disruption of the blood-brain barrier.\n   Higher GLCM contrast in MGMT-methylated cases may reflect\n   heterogeneous vascular architecture.\n\n3. SHAPE (ECCENTRICITY, SOLIDITY):\n   MGMT-unmethylated GBMs are often more infiltrative, leading to\n   irregular tumor shapes (lower solidity, higher eccentricity).\n   MGMT-methylated tumors may exhibit a more circumscribed morphology.\n\n4. SIGNAL INTENSITY (FLAIR MEAN, T2w MEAN):\n   MGMT methylation status influences the tumor's water content and\n   cellularity — reflected in T2/FLAIR signal levels. Higher FLAIR\n   mean often correlates with peritumoral edema extent.\n\n5. MULTI-MODALITY COMPLEMENTARITY:\n   - FLAIR: captures infiltrative peritumoral component\n   - T1wCE: reflects active angiogenesis and BBB breakdown\n   - T2w:   highlights edema and necrotic regions\n   - T1w:   baseline anatomy reference\n   Each provides a distinct biological view, hence fusion outperforms\n   any single modality.\n\n⚠️  STUDY LIMITATIONS & FUTURE DIRECTIONS\n──────────────────────────────────────────\n• N=50 is underpowered — full dataset (~585 patients) needed for\n  publication-quality conclusions\n• Pseudo-ROI (center crop + threshold) introduces noise —\n  use BraTS segmentation masks for true lesion-specific radiomics\n• Reproducibility: radiomics features depend on scan parameters;\n  harmonization (ComBat, normalization) needed for multi-site data\n• Deep learning approaches (CNNs, transformers) may outperform\n  hand-crafted features on this dataset\n• External validation cohort required before clinical translation\n\n📚 KEY REFERENCES\n──────────────────\n• Kickingereder et al. (2016) — Radiogenomics of GBM: MGMT & MRI\n• Aerts et al. (2014) — Decoding the tumor phenotype by noninvasive imaging\n• Chang et al. (2018) — Deep learning vs radiomics for MGMT prediction\n• Gillies et al. (2016) — Radiomics: images are more than pictures\n''')\n\nprint('═══════════════════════════════════════════════════════════')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:10:06.785112Z","iopub.execute_input":"2026-03-29T08:10:06.785469Z","iopub.status.idle":"2026-03-29T08:10:06.794111Z","shell.execute_reply.started":"2026-03-29T08:10:06.785438Z","shell.execute_reply":"2026-03-29T08:10:06.792631Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ─────────────────────────────────────────────────────────────\n# FINAL SUMMARY TABLE\n# ─────────────────────────────────────────────────────────────\n\nprint('\\n🏆 FINAL PERFORMANCE SUMMARY')\nprint('═' * 72)\nsummary = results_df.copy()\nsummary['AUC (mean±std)'] = summary.apply(\n    lambda r: f\"{r['AUC']:.3f} ± {r['AUC_std']:.3f}\", axis=1\n)\nprint(summary[['Features', 'Model', 'AUC (mean±std)', 'Accuracy', 'F1']]\n      .to_string(index=False, float_format='{:.3f}'.format))\n\nbest = results_df.loc[results_df['AUC'].idxmax()]\nprint(f'\\n🥇 Best configuration:')\nprint(f'   Model:    {best[\"Model\"]}')\nprint(f'   Features: {best[\"Features\"]}')\nprint(f'   AUC:      {best[\"AUC\"]:.3f} ± {best[\"AUC_std\"]:.3f}')\nprint(f'   Accuracy: {best[\"Accuracy\"]:.3f}')\n\nprint('\\n📁 Saved outputs:')\nfor fname in ['class_distribution.png', 'mri_modalities.png', 'preprocessing.png',\n              'roi_approximation.png', 'feature_selection.png', 'radiogenomic_association.png',\n              'modality_fusion.png', 'model_comparison.png', 'evaluation.png',\n              'interpretation_dashboard.png']:\n    print(f'   • {fname}')\n\nprint('\\n✅ Radiogenomics framework complete!')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-03-29T08:10:14.584539Z","iopub.execute_input":"2026-03-29T08:10:14.584898Z","iopub.status.idle":"2026-03-29T08:10:14.599995Z","shell.execute_reply.started":"2026-03-29T08:10:14.584866Z","shell.execute_reply":"2026-03-29T08:10:14.598742Z"}},"outputs":[],"execution_count":null}]}