TCGA-LIHC 肝癌转录组实战（上）

肝细胞癌（HCC，LIHC）是全球第六大常见肿瘤、第三大癌症死亡原因。TCGA-LIHC 队列一次性提供 373 例肿瘤 + 50 例癌旁正常的 HiSeqV2 转录组、配套的临床与生存数据，是验证『预后基因签名 + 分子分型』思路最经典的公开数据集之一。本系列两篇文章会把一份完整、可复现、可投稿的生信分析流程拆给你看 ——上篇从原始矩阵跑到一个判别力很强的 13 基因 LASSO-Cox 风险模型；下篇会把这个模型接到列线图、免疫微环境、WGCNA、分子亚型、通路活性、免疫治疗响应签名上，做完一整套机制 + 临床转化的故事线。

第一步 · 数据读入 + 样本注释 + PCA QC

我们直接读 UCSC Xena 镜像下来的三个文件：HiSeqV2 表达矩阵（log2 RSEM）、clinicalMatrix、Pan-Cancer survival。TCGA 的样本条码 14-15 位编码样本类型 —— 01-09 是肿瘤、10-19 是正常。表达矩阵做一次低表达过滤（要求 >25% 样本表达水平 > 1）后还剩 17 214 个基因 × 423 例样本。再用 MAD 取 top-2000 高变基因做 PCA。

R · 01_preprocess.R

# TCGA 条码 -> 样本类型
tcga_sample_type <- function(b) {
  code <- substr(b, 14, 15)
  ifelse(code %in% sprintf("%02d", 1:9), "Tumor",
  ifelse(code %in% sprintf("%02d", 10:19), "Normal", "Other"))
}
keep_g <- rowMeans(expr > 1) >= 0.25      # 17 214 基因
top_var <- names(sort(apply(expr_f, 1, mad), decreasing = TRUE))[1:2000]
pca <- prcomp(t(expr_f[top_var, ]), scale. = TRUE)

Top-2000 高方差基因 PCA — 肿瘤（红）与癌旁正常（蓝）在 PC1 上呈现非常干净的分离，PC1 解释方差约 34%，表明整体数据质量足够，无需额外的批次校正。

第二步 · limma 差异表达：3 653 个稳健 DEG

用 limma + trend + robust eBayes（log2 RSEM 数据直接走 trend 模式）做肿瘤 vs 癌旁正常对比。阈值取 |log2FC| > 1 且 FDR < 0.05，得到上调 1 605、下调 2 048 共 3 653 个差异基因。

R · 02_DEG.R

design <- model.matrix(~ 0 + samp$type)
colnames(design) <- levels(samp$type)            # Normal / Tumor
contr  <- makeContrasts(TvsN = Tumor - Normal, levels = design)
fit    <- lmFit(expr_f, design) |> contrasts.fit(contr) |>
          eBayes(trend = TRUE, robust = TRUE)
deg    <- topTable(fit, coef = "TvsN", number = Inf)
deg$dir <- with(deg, ifelse(adj.P.Val<0.05 & logFC> 1, "Up",
                     ifelse(adj.P.Val<0.05 & logFC<-1, "Down", "NS")))

Tumor vs Normal 差异基因火山图。红色为上调（n = 1 605）、蓝色为下调（n = 2 048）。标签为各方向 P 值 top 12 的代表基因。

Top 100 DEG（50 上 + 50 下）z-score 热图。样本按肿瘤 / 正常排序，可见两组之间几乎是反向的染色模式 —— 差异稳定且系统性。

第三步 · GO / KEGG / Hallmark 富集：肝代谢崩塌 + 细胞周期暴涨

GO 与 KEGG 用 clusterProfiler 的超几何检验分别对上调 / 下调集合做 ORA；Hallmark 通路则用全谱排序的 t 统计量跑 GSEA，既能看显著条目，也能看作用方向。

3.1 GO ORA · 三本体气泡图

GO 富集 top 5 / 本体 / 方向。下调集合（Down，蓝点）几乎完全被肝特异代谢与有机酸代谢通路占据；上调集合（Up，红点）则集中在染色体分离 / 有丝分裂 / DNA 复制等增殖通路。

3.2 KEGG ORA

KEGG 通路富集 top 10 / 方向。下调侧的 Drug metabolism / Retinol metabolism / Fatty acid degradation 直接对应了肝脏正常生理功能；上调侧的 Cell cycle / DNA replication 对应肿瘤增殖。

3.3 Hallmark GSEA · NES 排序

Hallmark GSEA 按 |NES| 排序的 top 20 通路。红色 = 在肿瘤中被激活、蓝色 = 在肿瘤中被抑制。E2F_TARGETS、G2M_CHECKPOINT、MYC_TARGETS、MITOTIC_SPINDLE 等增殖通路占据正向 top；BILE_ACID_METABOLISM、FATTY_ACID_METABOLISM、XENOBIOTIC_METABOLISM 等肝生理通路占据负向 top。

第四步 · 生存分析：从 5000 高变基因到 571 个 Cox 显著基因

在 367 例带 OS 信息的肿瘤样本上，对 top-5000 MAD 高变基因逐个跑 univariate Cox，BH 校正后得到 571 个 FDR < 0.05 的预后基因，与 DEG 取交集得到 432 个『又差异又预后』的稳健候选 —— 这就是下一步 LASSO 的输入特征池。

R · 04_survival.R

g_test <- names(sort(gene_mad, decreasing = TRUE))[1:5000]
uv <- do.call(rbind, lapply(g_test, function(g) {
  fit <- coxph(Surv(OS.time, OS) ~ as.numeric(xT[g, ]))
  s   <- summary(fit)
  data.frame(gene = g, HR = s$conf.int[,1],
             p = s$coefficients[, "Pr(>|z|)"])
}))
uv$fdr <- p.adjust(uv$p, method = "BH")            # 571 个 FDR<0.05

Univariate Cox top 20 基因森林图。红色 = 风险因子（HR>1），蓝色 = 保护因子（HR<1）。HR 用 log2 轴展示。

临床变量多因素 Cox 森林图。年龄、性别、grade 都不显著，只有 stage（III / IV）显著上调风险 —— 提示：在已有的临床变量里，可独立预测预后的只有分期。这正是我们后面要构建分子模型的动机。

按 AJCC stage 与组织学 grade 分层的 KM 曲线。Stage 之间分离明显（左），Grade 之间几乎重叠（右）—— 再次提示现有临床信息不够细。

第五步 · LASSO-Cox 风险模型：13 个基因，Train AUC > 0.80

把 432 个候选特征送进 10 折交叉验证的 LASSO-Cox，70 / 30 按事件分层切训练 / 测试，seed = 42 全程固定。λ.min 选出的 13 基因 panel：

R · 05_lasso_model.R

cvfit <- cv.glmnet(X_tr, y_tr, family = "cox",
                   alpha = 1, nfolds = 10,
                   type.measure = "deviance")
lam   <- cvfit$lambda.min
beta  <- coef(cvfit, s = lam)[, 1]
risk_score <- as.numeric(X %*% beta)
m$risk_group <- factor(ifelse(risk_score >= median_tr, "High", "Low"),
                       levels = c("Low", "High"))

5.1 LASSO 调参与系数路径

10-fold CV 偏差曲线（λ.min 与 λ.1se 用虚线标出）。λ.min 处选出 13 个非零系数。

LASSO 系数路径。λ 减小过程中各基因系数的进入次序，直观看出 IL18RAP / PSRC1 / SOCS2 等是最稳健的几个特征。

5.2 风险评分三联图

上：按风险评分排序的散点（红 = High，蓝 = Low）；中：每例患者的生存时间 + 终点状态（实心点为死亡）；下：top 20 LASSO 基因的 z-score 热图。三层叠加直观给出：高风险一侧死亡更密、肝代谢基因（蓝）更弱、增殖 / 应激基因（红）更强。

5.3 KM 与 time-ROC：训练 / 测试两个集合一致

Train 集 KM，High vs Low 风险组（log-rank p < 0.0001）。

Train 集 1/3/5 年 time-dependent ROC。AUC = 0.840 / 0.819 / 0.766。

Test 集 1/3/5 年 ROC。AUC = 0.738 / 0.738 / 0.735，三年时间窗维持稳定。

上篇小结 + 下篇预告

到这里上篇闭环：从原始 TCGA 矩阵 → 3 653 DEG → 571 个 Cox 显著基因 → 432 候选特征 → 13 基因 LASSO-Cox 风险模型。下篇会进一步：

① 把风险评分接进 5 变量列线图，做 1/3/5 年 OS 概率预测 + 校准曲线 + C-index 横向对比；

② ESTIMATE + 28 种免疫细胞 ssGSEA + 检查点表达，看高 / 低风险组的免疫微环境差异；

③ WGCNA 找出与风险评分相关的共表达模块和 hub 基因；

④ ConsensusClusterPlus 鉴定分子亚型 + 亚型 vs 临床 / 通路 / 免疫治疗签名比较。

如果你正打算用 TCGA 任意一个癌种重复这套思路，或者已经做了一半卡在某一步 —— 可以把数据丢给我们，复现 + 出图 + 报告一条龙。下篇见。

—— 长期承接 · 欢迎咨询 ——