406_放射组学Radiomics与基因组关联¶

# ========== 放射组学特征提取与基因组关联示例 ==========
import numpy as np
import pandas as pd
from scipy import stats
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

# ========== 1. 模拟放射组学特征 ==========
np.random.seed(42)
n_patients = 100  # 患者数
n_radiomic_features = 200  # 放射组学特征数（实际可达1000+）

# 生成放射组学特征矩阵（患者 x 特征）
radiomic_features = pd.DataFrame(
    np.random.randn(n_patients, n_radiomic_features),
    columns=[f"Radiomic_{i}" for i in range(n_radiomic_features)]
)

# ========== 2. 模拟基因组数据 ==========
n_genes = 500  # 基因数

# 基因表达数据（患者 x 基因）
gene_expression = pd.DataFrame(
    np.random.lognormal(3, 1, size=(n_patients, n_genes)),
    columns=[f"Gene_{i}" for i in range(n_genes)]
)

# 模拟EGFR突变状态（二分类）
egfr_status = pd.Series(
    np.random.binomial(1, 0.3, n_patients),  # 30%突变率
    name="EGFR_mutation"
)

# ========== 3. 放射组学特征提取（PyRadiomics示例框架） ==========
# 实际使用时需要安装：pip install pyradiomics
# import radiomics
# from radiomics import featureextractor

# 这里演示特征提取后的处理流程
# 真实代码：
# extractor = radiomics.featureextractor.RadiomicsFeatureExtractor()
# result = extractor.execute(image_path, mask_path)  # 输入影像和分割掩码

# ========== 4. 特征选择（去除高度相关特征） ==========
# 计算特征间相关系数矩阵
corr_matrix = radiomic_features.corr().abs()  # 取绝对值

# 找出相关系数>0.9的特征对（高度冗余）
upper_tri = corr_matrix.where(
    np.triu(np.ones(corr_matrix.shape), k=1).astype(bool)  # 只看上三角
)
highly_correlated = [
    col for col in upper_tri.columns 
    if any(upper_tri[col] > 0.9)  # 与任何特征相关系数>0.9
]
print(f"高度相关特征数（待删除）: {len(highly_correlated)}")

# 删除冗余特征
radiomic_filtered = radiomic_features.drop(columns=highly_correlated)
print(f"过滤后特征数: {radiomic_filtered.shape[1]}")

# ========== 5. 放射组学特征与EGFR突变关联 ==========
# 对每个放射组学特征做Mann-Whitney U检验（比较突变vs野生型）
association_results = []

for feature in radiomic_filtered.columns:
    mutant_values = radiomic_filtered[feature][egfr_status == 1]    # EGFR突变患者
    wildtype_values = radiomic_filtered[feature][egfr_status == 0]  # 野生型患者

    # Mann-Whitney U检验（不假设正态分布）
    u_stat, p_value = stats.mannwhitneyu(
        mutant_values, wildtype_values, alternative="two-sided"
    )

    # 计算效应量（AUC = U/(n1*n2)，AUC>0.7表示较强区分能力）
    auc = u_stat / (len(mutant_values) * len(wildtype_values))
    auc = max(auc, 1 - auc)  # 确保AUC>0.5

    association_results.append({
        "feature": feature,
        "p_value": p_value,
        "AUC": auc
    })

assoc_df = pd.DataFrame(association_results)

# FDR校正
from scipy.stats import rankdata
p_vals = assoc_df["p_value"].values
ranks = rankdata(p_vals)
fdr = np.minimum(p_vals * len(p_vals) / ranks, 1.0)
assoc_df["FDR"] = fdr

# 筛选显著关联特征
sig_features = assoc_df[(assoc_df["FDR"] < 0.05) & (assoc_df["AUC"] > 0.65)]
print(f"\n与EGFR突变显著关联的放射组学特征数: {len(sig_features)}")

# ========== 6. 构建预测模型 ==========
from sklearn.linear_model import LogisticRegressionCV
from sklearn.model_selection import cross_val_score

# 用显著特征构建EGFR突变预测模型
if len(sig_features) > 0:
    X = radiomic_filtered[sig_features["feature"].tolist()].values
    y = egfr_status.values

    # 标准化
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)

    # Lasso逻辑回归（自动特征选择）
    model = LogisticRegressionCV(
        cv=5,            # 5折交叉验证
        penalty="l1",   # Lasso正则化，自动选特征
        solver="liblinear",
        random_state=42
    )

    # 5折交叉验证评估AUC
    cv_scores = cross_val_score(model, X_scaled, y, cv=5, scoring="roc_auc")
    print(f"预测EGFR突变的AUC: {cv_scores.mean():.3f} ± {cv_scores.std():.3f}")