TCGA-LUAD 多组学预后建模实战（上）

肺腺癌（Lung Adenocarcinoma，LUAD）是非小细胞肺癌中发病最高的亚型，尽管 EGFR-TKI 与免疫检查点抑制剂已显著改善部分患者预后，但相当一部分人在两年内仍会复发或转移。围绕 LUAD 建立一套可解释、可复现的预后基因模型，对于早期分层、术后随访以及辅助治疗决策都具有非常直接的临床价值。本文是 LUAD 系列的上篇，使用 TCGA-LUAD 的 HiSeqV2 mRNA、miRNA、临床矩阵与生存数据，依次走完 5 个标准步骤：① 数据预处理；② 肿瘤 vs 正常差异分析；③ 单因素 Cox → LASSO-Cox → 风险评分 → 时间依赖 ROC → 列线图；④ 风险评分与临床变量的关联；⑤ GO / KEGG / GSEA 通路富集。

第一步 · 数据预处理：四类原始文件合成一个分析对象

使用 data.table::fread 一次性读入 mRNA 表达矩阵、miRNA 表达矩阵、临床矩阵与生存表，通过 TCGA Barcode 第 14-15 位标注 Tumor / Normal，再以 patient ID 把生存 / 临床 / 表达对齐到同一份样本表。预处理产物 luad_data.rds 是后续所有脚本的统一输入。

R · 01_data_preprocessing.R

mrna <- fread(f_mrna,  data.table = FALSE)
mir  <- fread(f_mir,   data.table = FALSE)
clin <- fread(f_clin,  data.table = FALSE,
              na.strings = c("", "NA", "[Not Available]", "[Unknown]"))
surv <- fread(f_surv,  data.table = FALSE)

mrna_type <- tcga_sample_type(colnames(mrna))   # 14-15 位 -> Tumor / Normal
tumor_ids <- intersect(colnames(mrna)[mrna_type == "Tumor"], surv$sample)
pheno     <- build_pheno(tumor_ids, surv, clin)

1.1 队列总览

Figure 1. TCGA-LUAD 队列概览。从左上至右下分别为：mRNA / miRNA 各组别样本数；病理分期分布（I–IV）；诊断年龄分布；以及随访时长与生存状态。样本总量、肿瘤/正常配对、Stage 与年龄分布均符合后续建模的基本要求。

第二步 · 差异表达：肿瘤 vs 正常（mRNA + miRNA）

差异分析使用 limma + trend。阈值 |log2FC| > 1 且 FDR < 0.05，分别对 mRNA 与 miRNA 做 Tumor vs Normal。

R · 02_differential_expression.R

design <- model.matrix(~ grp)   # grp: Normal, Tumor
fit    <- eBayes(lmFit(expr, design), trend = TRUE)
tt     <- topTable(fit, coef = 2, number = Inf, sort.by = "P")
tt$Direction <- with(tt,
   ifelse(adj.P.Val < 0.05 & logFC >  1, "Up",
   ifelse(adj.P.Val < 0.05 & logFC < -1, "Down", "NS")))

两个组学的差异分布如下（数量级与典型 TCGA 报道一致）：

Figure 2. mRNA / miRNA 差异基因汇总。mRNA 共筛得 1820 个上调 + 2610 个下调；miRNA 102 上调 + 32 下调。上调以增殖、糖代谢、细胞周期基因为主，下调以肺组织特异分化、免疫相关基因为主。

2.1 mRNA 火山图（差异最显著的基因被标注）

Figure 3. mRNA volcano plot. 右上深色点为上调 DEG，左上为下调 DEG。标注的是 FDR 最显著的前 15 个 DEG，包括 SPP1、CRABP2、KRT6A、MAGE 家族等。

2.2 miRNA 火山图

Figure 4. miRNA volcano plot. miR-21、miR-210、miR-183 等经典 oncomiR 显著上调；let-7 家族（let-7a / let-7c）作为典型抑癌 miRNA 显著下调 —— 与文献完全一致。

2.3 Top DEG 热图

Figure 5. Top 50 上调 + 50 下调 mRNA 热图。Tumor 与 Normal 在 Z-score 层面形成清晰的对比块，说明 DEG 集合具有强分类能力。

第三步 · 预后建模：单因素 Cox → LASSO → 风险评分

建模流程严格遵循『先筛选、再压缩、再外验』的三段式：首先对每个候选基因做单因素 Cox，保留 p < 0.01 的稳健预后基因；再用 10 折交叉验证的 LASSO-Cox 把多重共线性压成稀疏权重；用所选基因的线性组合计算每位患者的连续 risk score；以中位数划高低风险组，做 KM、时间依赖 ROC 与列线图。

R · 03_survival_analysis.R

# 单因素 Cox 筛预后基因
uni <- lapply(candidate, function(g) coxph(Surv(OS.time, OS) ~ expr[g, ],
                                            data = pheno))
# LASSO-Cox + CV 选 lambda.min
cvfit <- cv.glmnet(X, Surv(OS.time, OS), family = "cox", alpha = 1, nfolds = 10)
sel   <- coef(cvfit, s = "lambda.min")              # 稀疏权重
risk  <- as.numeric(X %*% sel)                      # 连续 risk score

3.1 顶部预后基因（单因素 Cox Top 20）

Figure 6. 单因素 Cox 前 20 名预后基因 Forest plot。右侧（HR > 1，coral）为风险基因，左侧（HR < 1，teal）为保护基因。

3.2 LASSO 系数路径与最优 lambda

LASSO 系数路径 + CV deviance — Figure 7. 左：LASSO 系数路径，每一条曲线代表一个基因的权重随 log(λ) 的衰减；右：10 折交叉验证的 partial likelihood deviance，coral 虚线为 lambda.min。

3.3 风险评分与生存状态全景

Risk score 排序 — Figure 8. Risk score 按升序排列，每个点为一位患者。Coral = 高风险组（risk > 中位数），Teal = 低风险组。

生存时间 + 生死状态 — Figure 9. 与 Figure 8 同序，纵轴为生存时间，颜色为生死状态。高风险段（右侧）死亡事件密度明显上升，与 risk score 升序方向一致。

13 个 LASSO 基因 Z-score 热图 — Figure 10. 13 个 LASSO 基因在样本顺序上的 Z-score 热图。高风险组多数 LASSO 基因表达整体上调，可视化地解释了 risk score 的来源。

3.4 KM 生存曲线：高低风险组

Figure 11. 风险分组的 Kaplan-Meier 曲线。高风险组生存显著差于低风险组（log-rank p << 0.001），分组效应非常稳定。

3.5 关键单基因 KM

Figure 12. RGS20 单基因 KM。高表达组生存显著较差，与 LASSO 系数为正方向一致。

3.6 时间依赖 ROC（1 / 3 / 5 年）

3.7 Nomogram + 3 年校准曲线

Figure 16. 联合 risk score + Stage + 年龄的 Nomogram，可对单个患者直读 1 / 3 / 5 年死亡概率。

第四步 · 风险评分与临床变量的关联

一个好的预后评分应该和肿瘤分期、淋巴结、转移这些既有临床指标在方向上保持一致；同时不应被性别等无关变量主导。我们对每个临床变量做 Kruskal-Wallis 检验，只保留显著的图。

4.1 风险评分随病理分期升高

Figure 18. Risk score by Pathologic Stage (I-IV). Kruskal-Wallis p < 1e-10. Stage 越高，risk score 越高，方向完全符合临床直觉。

4.2 T / N / M 分期

Figure 19. Risk score by Pathologic T. T 分期越高 risk score 越高 (p ≈ 1e-5)。

Figure 20. Risk score by Pathologic N. 淋巴结转移阳性者 risk score 显著升高 (p ≈ 1e-7)。

Figure 21. Risk score by Pathologic M. M1 患者 risk score 较 M0 显著升高 (p ≈ 4e-3)。

4.3 吸烟史

Figure 22. Risk score by tobacco smoking history. 重度吸烟史 (4/5) 组 risk score 系统性偏高 (p ≈ 5e-3)，与流行病学一致。

4.4 LASSO 顶级基因在不同 Stage 间的表达

Top 6 LASSO 基因在 Stage I-IV 表达分布 — Figure 23. Top 6 LASSO 基因在 Stage I-IV 间的表达分布。多数风险基因随 Stage 升高单调上升，模型与肿瘤分期生物学一致。

4.5 单因素 + 多因素 Cox：Risk score 是独立预后因子

Figure 24. Univariate + Multivariate Cox Forest plot。Risk score 在调整 age / Stage / gender 之后仍为独立预后因子，HR 显著大于 1。

第五步 · 功能富集：GO + KEGG + GSEA

用 clusterProfiler 对上 / 下调 DEG 分别做超几何检验（ORA），再以 limma 的 t 统计量作 ranked list 跑 GSEA，两种视角互为印证。

R · 05_enrichment_analysis.R

# ORA
go <- enrichGO(gene = ent, OrgDb = org.Hs.eg.db, ont = "BP",
               pvalueCutoff = 0.05, readable = TRUE)
kegg <- enrichKEGG(gene = ent, organism = "hsa", pvalueCutoff = 0.05)

# GSEA on ranked t-statistic
glist <- setNames(sort(deg$t, decreasing = TRUE), deg$ENTREZ)
gsea_kegg <- gseKEGG(geneList = glist, organism = "hsa",
                     pvalueCutoff = 0.05)

5.1 GO-BP：上调 DEG 富集

Figure 25. GO-BP enrichment of up-regulated DEGs. 前列条目以细胞周期、DNA 复制、纺锤体组装、染色体分离为主，提示肿瘤组织高度增殖的表型。

5.2 GO-BP：下调 DEG 富集

Figure 26. GO-BP enrichment of down-regulated DEGs. 肺组织发育、纤毛运动、表面活性物质合成、肺泡分化 —— 全部是肺正常生理功能的丢失。

5.3 KEGG 通路富集（上 / 下调）

Figure 27. KEGG enrichment of up-regulated DEGs. 细胞周期、p53 信号、DNA 复制位列前茅。

Figure 28. KEGG enrichment of down-regulated DEGs. 钙信号、cAMP 信号、神经活性配体受体相互作用等正常组织通路富集。

5.4 GSEA · GO-BP 双向气泡图

Figure 29. GSEA on GO BP. 左侧为正向（在肿瘤中上调）富集通路，右侧为负向。两侧通路与 ORA 结论一致且更细致。

5.5 GSEA · KEGG ridge

Figure 30. GSEA on KEGG, ridge plot. 富集通路按 enrichment score 排序，每条曲线为该通路 leading edge 基因的 logFC 密度。

5.6 GSEA · 单通路 enrichment plot

最显著的 KEGG 通路（细胞周期）GSEA — Figure 31. 最显著的 KEGG 通路（细胞周期）的 GSEA enrichment plot。

Figure 32. 第二个显著 KEGG 通路的 enrichment plot。

上篇小结 · 下篇预告

上篇我们已经把 LUAD 的『从矩阵到风险评分到通路』这条主线走完了：504 例配对样本、4430 个 DEG、13 基因 LASSO 模型、Stage 关联、GO/KEGG 富集 —— 一个完整、可投稿的预后分析骨架。

下一篇会把 WGCNA 共表达模块 / hub 基因、ESTIMATE + ssGSEA 免疫浸润 + 22 个免疫检查点、Hallmark 50 通路单样本 GSVA 评分、分子亚型（ConsensusClusterPlus）、C-index 自助法 + DCA + 随机生存森林反向验证、miRNA-mRNA 调控网络（multiMiR）、以及驱动基因 OncoPrint + Stage→Cluster→Risk→Status 流向图，一次性补齐 LUAD 论文的所有进阶分析模块。