肿瘤预后模型构建与验证¶

特征筛选 → 模型构建 → 内部验证 → 外部验证 → 临床应用

# === LASSO-Cox 预后模型完整流程 ===
library(glmnet)                        # LASSO/弹性网络
library(survival)                      # 生存分析
library(survminer)                     # 生存曲线可视化
library(timeROC)                       # 时间依赖ROC

# 1. 准备数据
expr <- read.table("expression_tpm.txt",  # 表达矩阵
                   header = TRUE, row.names = 1)
clinical <- read.csv("clinical.csv")   # 临床信息
# 确保样本匹配
common <- intersect(colnames(expr), clinical$sample_id)
expr <- expr[, common]                 # 子集表达
clinical <- clinical[clinical$sample_id %in% common, ]

# 2. 特征预筛选（单因素Cox过滤）
genes <- rownames(expr)                # 所有基因
uni_results <- data.frame(             # 单因素结果
  gene = character(), pvalue = numeric(), HR = numeric()
)

for (gene in genes) {                  # 遍历每个基因
  gene_expr <- as.numeric(expr[gene, ])  # 该基因表达量
  surv_obj <- Surv(clinical$OS_time, clinical$OS_status)  # 生存对象
  fit <- coxph(surv_obj ~ gene_expr)   # 单因素Cox
  p <- summary(fit)$coefficients[, 5]  # p值
  hr <- summary(fit)$coefficients[, 2] # HR
  uni_results <- rbind(uni_results,
    data.frame(gene = gene, pvalue = p, HR = hr))
}

# 筛选 p < 0.05 的基因进入 LASSO
sig_genes <- uni_results$gene[uni_results$pvalue < 0.05]
cat("单因素显著基因数:", length(sig_genes), "\n")

# 3. LASSO-Cox 变量选择
x <- t(expr[sig_genes, ])              # 转置为 样本×基因
y <- Surv(clinical$OS_time, clinical$OS_status)  # 生存数据

# 10折交叉验证选择最优lambda
set.seed(42)                           # 固定随机种子
cv_fit <- cv.glmnet(                   # 交叉验证
  x = as.matrix(x), y = y,            # 输入
  family = "cox",                      # Cox模型
  alpha = 1,                           # LASSO (alpha=1)
  nfolds = 10                          # 10折
)

plot(cv_fit)                           # 绘制lambda-偏差曲线
best_lambda <- cv_fit$lambda.min       # 最优lambda

# 提取选中的基因
coef_matrix <- coef(cv_fit, s = best_lambda)  # 系数
selected <- rownames(coef_matrix)[coef_matrix[, 1] != 0]
cat("LASSO选中基因数:", length(selected), "\n")
cat("基因:", paste(selected, collapse = ", "), "\n")

# 4. 构建最终模型（多因素Cox）
model_data <- data.frame(              # 建模数据
  time = clinical$OS_time,             # 生存时间
  status = clinical$OS_status,         # 生存状态
  t(expr[selected, ])                  # 选中基因表达
)
final_model <- coxph(                  # 多因素Cox
  Surv(time, status) ~ .,             # 所有变量
  data = model_data
)
summary(final_model)                   # 查看模型

# 5. 计算风险评分
risk_score <- predict(final_model,     # 预测
                      type = "lp")     # 线性预测值
model_data$risk <- ifelse(             # 分组
  risk_score > median(risk_score),     # 高于中位数
  "High", "Low")

# 6. Kaplan-Meier 可视化
km_fit <- survfit(Surv(time, status) ~ risk,
                  data = model_data)   # KM拟合
ggsurvplot(km_fit,                     # KM曲线
           pval = TRUE,               # 显示p值
           risk.table = TRUE,          # 风险表
           palette = c("red", "blue")) # 颜色

# 7. 时间依赖 ROC
roc <- timeROC(T = model_data$time,    # 生存时间
               delta = model_data$status,  # 事件状态
               marker = risk_score,    # 风险评分
               times = c(365, 1095, 1825),  # 1/3/5年
               cause = 1)
plot(roc, time = 1095)                 # 3年ROC曲线
cat("1年AUC:", roc$AUC[1], "\n")
cat("3年AUC:", roc$AUC[2], "\n")
cat("5年AUC:", roc$AUC[3], "\n")