跳到正文
百沐一下百沐一下
登录
本节目录

复现 TCGA 泛癌预后分析

跟完这篇，您能挑一篇「某基因在 TCGA 里泛癌预后」的论文，把它的取数、按癌种分组、逐癌种生存分析和那张跨癌种森林图完整走一遍，最后得到一张和原文趋势一致的森林图——并且在对不上的时候，知道该去哪里找原因。

泛癌预后图是肿瘤生信里最高频的一类图：一个基因，横跨几十个癌种，一图看尽它在每个癌里「表达高低和病人活得久不久」的关系。它套路固定、数据公开，是练「复现」再合适不过的靶子。这篇就带您把这套路完整走一遍。

先抄准论文的分析口径

动手之前，先把论文的方法部分（Methods）和图注读透。泛癌预后能不能对上，几乎全取决于下面这几项抄得准不准。建议边读边填一张表：

要素在论文哪里找为什么关键
目标基因标题 / 摘要 / Methods复现的对象
纳入哪些癌种Methods、图是 33 个全癌，还是剔掉了样本 / 事件太少的子集
数据来源Data availabilityTCGA / GDC 还是 UCSC Xena，决定取数路径
表达单位Methodslog2(TPM+1) / FPKM / counts，影响分组切点
生存终点Methods、图注OS / PFI / DSS / DFI，换个终点结论可能就变
分组方式Methods按中位数二分？最佳截断？三分位？
生存方法MethodsKM + log-rank，还是单因素 Cox 出 HR
森林图指标图注一般是 HR + 95% 置信区间
多重校正Methods跨癌种有没有做 BH 校正、阈值定在哪

几个读论文时容易漏、却很要命的细节：

  1. 生存终点到底用的哪个。 TCGA 常见四种终点：总生存 OS、无进展间期 PFI、疾病特异生存 DSS、无病间期 DFI。它们的 time 和 event 定义不同，同一个基因换个终点，HR 和显著性都可能变。论文用 OS，您就跟 OS，别混。
  2. 分组是中位数还是「最佳截断」。 按中位数二分最稳、最好复现；用最佳截断（surv_cutpoint）是从结局里挑切点，天然让 p 偏小，换个数据就难重现。论文用哪种，您就用哪种。
  3. 纳入了哪些癌种、剔了哪些样本。 有的论文用全 33 癌，有的会去掉事件数太少的癌种，或只留原发肿瘤样本。漏掉这句，您每个癌种的样本数就和人家对不上。

复现的目标是趋势一致、结论可重现，不是像素级一模一样。只要森林图里「哪些癌种是危险因素（HR>1）、哪些是保护因素」的大格局对得上，主要癌种的 HR 量级也接近，就算复现成功。追求每个 HR 分毫不差，多数时候既做不到、也没必要。

取到 TCGA 表达与临床终点

泛癌预后要两样东西：一张跨癌种的表达矩阵，一张带生存终点的临床表。您几乎不用从原始测序数据算起，主流有两个入口。

  • UCSC Xena（推荐上手）：UCSC Xena 把 TCGA 泛癌（PanCanAtlas）的表达整理成了一张规整、已归一化、开箱即用的大矩阵（值通常是 log2(norm_count+1) 这类），并配套提供一份整理好的生存 / 表型表，里面就有 OS、PFI、DSS、DFI 及各自的时间列。取数最省事，优先走这条。
  • GDC / TCGAbiolinks：GDC 是 TCGA 的官方门户，字段最全、最原始，但整理和统一的活得自己干。需要最原始数据时再回到这里。

生存终点表的列名通常长这样，.time 是随访时长（天）、不带后缀的是 0/1 事件：

列含义惯例建议
OS / OS.time总生存：是否死亡 / 随访天数大多数癌种都相对可靠，最常用
PFI / PFI.time无进展间期随访较短的癌种里比 OS 更灵敏，也常用
DSS / DSS.time疾病特异生存依赖死因标注，部分癌种数据质量有限
DFI / DFI.time无病间期可用样本较少，用前先看覆盖度

（OS 与 PFI 通常是最稳妥的两个终点，DSS、DFI 在部分癌种里数据不全，用哪个以论文为准。）

# 从 UCSC Xena 下载两份文件(具体文件名以 Xena 页面为准):
# 1) TCGA 泛癌表达矩阵, 值多为 log2(norm_count+1), 行基因、列样本
# 2) 配套生存/表型表, 含 OS、OS.time、PFI、PFI.time、cancer type 等列
head TCGA_PanCan_expression.tsv     # 先看一眼形态: 行是基因、列是样本?
head TCGA_survival.tsv              # 终点列名、时间单位对上没有?

取回后在 R 里把两张表按样本 ID 对齐，并按论文口径过滤样本——最常见的是只留原发肿瘤（TCGA 样本条码第 4 段为 01）、去掉重复样本：

# expr: 表达矩阵(行基因、列样本); surv: 生存/表型表, 每行一个样本
gene <- "TP53"                                   # 换成您要复现的目标基因

dat <- data.frame(
  sample = colnames(expr),
  expr   = as.numeric(expr[gene, ]),
  stringsAsFactors = FALSE
)
dat <- merge(dat, surv, by = "sample")           # 按样本 ID 对齐表达与生存

# 只留原发肿瘤样本(条码第 4 段 01), 去重
dat <- dat[substr(dat$sample, 14, 15) == "01", ]
dat <- dat[!duplicated(dat$sample), ]

一个务实的提醒：能用 Xena 整理好的矩阵和终点表，就别自己从 GDC 原始数据拼。 泛癌复现要对齐几十个癌种的样本筛选、表达归一化、终点定义一大串口径，用现成的规整数据，您复现的是「从矩阵到森林图」这一段，反而更干净、更可控。数据下下来怎么整理成规范矩阵，见 表达矩阵整理。

逐癌种：分组 + 单因素 Cox

数据就位后，核心就一件事：在每个癌种里各跑一次生存分析，把结果攒成一张表。泛癌森林图的每一行，对应一个癌种的单因素 Cox 结果（HR、置信区间、p）。

按论文的分组方式来。若论文按表达中位数把每个癌种的病人切成高、低两组，就照做——注意中位数要在每个癌种内部各自算，不能全泛癌用一个切点：

library(survival)

cancers <- unique(dat$cancer)          # 各癌种缩写, 如 BRCA、LUAD、LIHC……

res <- lapply(cancers, function(ct) {
  d <- subset(dat, cancer == ct)
  d <- d[!is.na(d$OS.time) & !is.na(d$OS), ]     # 该终点缺失的样本先剔掉

  # 在该癌种内部按中位数二分, 以低表达为参照——参照方向全癌统一, 否则 HR 会取到倒数
  d$group <- ifelse(d$expr > median(d$expr, na.rm = TRUE), "high", "low")
  d$group <- factor(d$group, levels = c("low", "high"))

  fit <- coxph(Surv(OS.time, OS) ~ group, data = d)   # 论文用哪个终点就换哪个
  s   <- summary(fit)
  data.frame(
    cancer = ct,
    n      = nrow(d),
    HR     = s$conf.int[, "exp(coef)"],
    lower  = s$conf.int[, "lower .95"],
    upper  = s$conf.int[, "upper .95"],
    p      = s$coefficients[, "Pr(>|z|)"]
  )
})
res <- do.call(rbind, res)

读这张结果表：HR>1 表示在该癌种里目标基因高表达是危险因素（关联更差的生存），HR<1 是保护因素，置信区间不跨过 1 与 p<0.05 是一回事，都提示关联显著。

两个常见变体，按论文口径二选一：

  • 把表达当连续变量直接进 Cox（不二分）：coxph(Surv(OS.time, OS) ~ expr, data = d)。此时 HR 是「表达每升高一个单位」的效应，解释和二分不同，别混着报。
  • KM + log-rank 出 p、Cox 出 HR：想同时给生存曲线，就对代表性癌种用 survfit + ggsurvplot 画 KM，森林图仍用 Cox 的 HR。KM 与 log-rank 怎么读，见 生存分析入门；Cox 的 HR、参照水平、比例风险假设怎么讲究，见 Cox 回归。

阈值 p<0.05 是文献里最常见的默认，属于常用惯例、非铁律。但复现语境下您不用纠结哪套「更好」——照抄论文那套就对了。判断阈值本身合不合理，是原作者的事，不是复现要解决的问题。

汇成一张跨癌种森林图

每个癌种一行结果攒齐，就能拼成泛癌森林图：一行一个癌种，点是 HR，横线是 95% 置信区间，竖直参考线在 HR=1——线整条落在 1 右侧且不压到 1，说明在该癌里是确凿的危险因素；压过 1 则不显著。

library(ggplot2)

# 按 HR 排序, 图看起来更有层次
res$cancer <- factor(res$cancer, levels = res$cancer[order(res$HR)])

ggplot(res, aes(x = HR, y = cancer)) +
  geom_point(size = 2) +
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.25) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "grey40") +
  scale_x_log10() +          # HR 惯例用对数轴, 让 0.5 与 2 在视觉上对称
  labs(x = "Hazard ratio (95% CI)", y = NULL) +
  theme_bw()

对着原图逐项校准：

  • 癌种排序：原图按 HR 排就按 HR 排，按字母序就按字母序，好逐行对照。
  • 对数横轴：HR 用 scale_x_log10() 是通行做法，保护因素（HR<1）和危险因素（HR>1）在图上对称，别用线性轴把 HR=0.5 压成一小段。
  • 显著性标注：原图常用星号或加粗标出显著癌种，把您的 p（或下一节的 padj）对应标上即可。

想要更强的排版控制（HR 与 CI 数值列排在图旁、分组标题），可用 forestplot 包手动搭；快速出图用上面的 ggplot2 足够。

校正多重比较，别被随机撞出的显著骗了

这一节很容易被跳过，却是泛癌分析的硬伤高发区。同一个基因在 33 个癌种里各做一次检验，纯靠运气也会撞出几个 p<0.05。 直接数「有几个癌种显著」而不校正，几乎必然高估。

标准做法是对所有癌种的 p 值做一次多重检验校正（BH 法最常用），看校正后还剩几个站得住：

res$padj <- p.adjust(res$p, method = "BH")   # 跨癌种一起校正
sum(res$padj < 0.05)                          # 校正后仍显著的癌种数

几点要留心：

  • padj<0.05 是常用惯例、非铁律，按检验数量和研究目的调。复现时以论文的口径为准：论文若做了 BH，您也做；论文报的是原始 p，您对齐原始 p，但心里清楚它偏乐观。
  • 小样本癌种的结论天然更不稳。 某些癌种只有几十例、事件数个位数，Cox 会给出很宽的置信区间、飘忽的 HR。森林图里这类行「点看着极端、线宽得吓人」，解读要留足余地。
  • 相关不是因果。 森林图给的是「表达高低与生存的关联」，不是「这个基因驱动了预后」。措辞老实地停在「关联」「提示」。

对不上怎么排查

第一次复现几乎不会一次对齐——这很正常。按下面几类原因从大到小逐个排除，它们对结果的影响通常也是从上往下递减：

可能原因典型表现怎么查、怎么办
生存终点选错整张森林图方向或显著格局都对不上核对论文用的是 OS 还是 PFI / DSS / DFI，换成同一个终点
参照方向反了所有 HR 恰好取成了倒数（0.5↔2）统一以低表达为参照，或统一「高 vs 低」的方向
样本筛选不同每个癌种的样本数和论文对不上查是否只留原发肿瘤（01）、去重复、按论文剔除了哪些癌种
分组切点不同部分癌种显著性翻转对齐「每癌内部按中位数」还是最佳截断；别全泛癌用一个切点
表达单位不同HR 量级系统性偏移对齐 log2(TPM+1) 还是 counts / FPKM
数据版本不同结果整体轻微漂移对齐 Xena / GDC 的数据版本快照，记下下载日期
多重校正口径「显著癌种数」差很多复查论文是 BH 校正后还是原始 p、阈值是 padj 还是 p

两条实用心法：

  1. 从影响最大的往下查。 终点、参照方向、样本筛选这类「结构性」差异会让整张图面目全非；数据小版本差异只带来毫厘漂移。先排大的，再抠小的。
  2. 记录您的环境。 跑完 sessionInfo() 把 R 和 survival 等包的版本存下来，既方便自己回溯，也是别人将来复现您结果的依据。

诚实地说：想和几年前的论文分毫不差，往往做不到——当年的 TCGA 数据快照、样本筛选口径可能已经很难完整还原。这不是您做错了。复现追求的是「危险 / 保护的大格局和主要癌种的结论能重现」，把差异查清、能解释清楚，本身就是有价值的产出。

在平台上跑一遍

上面这套「取 TCGA 数据 → 逐癌种分组 → 单因素 Cox → 跨癌种森林图 → 多重校正」，在 百沐一下 不写代码也能走完：说清目标基因、生存终点和分组方式（比如「对 TP53 做泛癌预后分析，用 OS 终点、每个癌种按表达中位数分高低组，逐癌种单因素 Cox，画森林图并做 BH 校正」），平台会从公共数据库取数、分癌种跑完生存分析、画出森林图，并连同可复现的 R 代码与用到的包版本一起交付，方便您和原文逐项比对。想连自己的数据一起分析也可以，直接上传即可。完整流程见 组学分析助手使用指南。

相关