TCGA数据的单基因预后初筛实战：Cox+KM

一、写在前面

面向生信初学者的一篇实操合集，分成两个独立部分：第一部分教你用TCGAbiolinks从GDC下载TCGA表达数据；第二部分教你拿到数据后，对你自己选定的基因做单因素Cox回归 + 中位数截断KM生存分析（以FSTL3为例，把GENE <- "FSTL3"换成自己的基因名即可）。

需要说明的是，第二部分完成的实际上是在 TCGA 队列中的单基因预后关联分析（初筛），而不是严格意义上的”独立外部验证”。

二、TCGA 数据下载

TCGA数据获取最常用的方案是TCGAbiolinks，核心四个函数：查询 → 下载 → 整理 → 保存。

2.1 查询：GDCquery

library(TCGAbiolinks)library(SummarizedExperiment)query<-GDCquery(project="TCGA-LUAD",# 项目 IDdata.category="Transcriptome Profiling",# 转录组data.type="Gene Expression Quantification",# 基因表达定量workflow.type="STAR - Counts",# STAR 流程（TCGA 官方推荐）sample.type="Primary Tumor"# 只要原发肿瘤)

运行输出如下（能打印出o Preparing output即表示查询成功）：

-------------------------------------- o GDCquery: SearchinginGDC database -------------------------------------- Genome of reference: hg38 -------------------------------------------- oo Accessing GDC. This might take a while... -------------------------------------------- ooo Project: TCGA-LUAD -------------------- oo Filtering results -------------------- ooo By data.type ooo By workflow.type ooo By sample.type ---------------- oo Checking data ---------------- ooo Checkingifthere are duplicated cases ooo Checkingifthere are resultsforthe query ------------------- o Preparing output -------------------

📌 通俗理解：GDCquery 就像在淘宝里按条件搜索商品——“我要 LUAD 的转录组定量、STAR 流程、原发肿瘤样本”，它返回一张”购物清单”，而不是数据本身。

🚫 务必避坑：sample.type别漏了。如果省略 sample.type，会同时下载 Solid Tissue Normal（癌旁正常组织）和 Primary Tumor 等所有类型。如果研究目标是肿瘤患者预后分析，建议在查询阶段明确限定 Primary Tumor，避免把不同样本类型混入预后队列。

2.2 下载：GDCdownload

# 下载（api 方式，每次 20 个文件，断点续传更稳） GDCdownload(query, method = "api", files.per.chunk = 20)

运行输出如下（会显示总文件数、总大小，并按 chunk 分批下载）：

Downloading dataforproject TCGA-LUAD GDCdownload will download540files. A total of2.288476735GB Downloading chunk1of27(20files, size=84.718875MB)as Sun_Sep__6_14_18_37_2026_0.tar.gz

2.3 查看下载的文件结构：list.files

下载完成后，可以用 list.files() 查看文件结构：

list.files("GDCdata",recursive=TRUE)# 递归列出所有文件length(list.files("GDCdata",recursive=TRUE))# 文件总数

list.files的结果如下（实际会列出全部 540 个文件，此处为目录结构示意）：

GDCdata/ └── TCGA-LUAD/ └── harmonized/ └── Transcriptome_Profiling/ └── Gene_Expression_Quantification/ ├── Sun_Sep__6_14_18_37_2026_0.tar.gz# 分块包├── Sun_Sep__6_14_18_37_2026_1.tar.gz └──...# 共 27 个 chunk

📌 提示：这些.tar.gz 是分块打包的原始文件。通常不需要手动逐个处理，后续可通过GDCprepare()对已下载数据进行整理、合并成 SummarizedExperiment对象。

2.4 整理与保存：GDCprepare + saveRDS

# 整理成 SummarizedExperiment se <- GDCprepare(query, summarizedExperiment = TRUE) # 保存为 R 原生对象（后续分析直接 readRDS） saveRDS(se, "TCGA-LUAD_se.rds")

GDCprepare得到的se是标准的SummarizedExperiment对象：assay() 存表达矩阵，colData() 存临床信息，rowData() 存基因注释。

💡 关键知识点：STAR - Counts 流程会给出三种定量方式，都在 assayNames(se)里：

assayNames(se)# [1] "unstranded" # 原始 counts# [2] "tpm_unstrand" # TPM# [3] "fpkm_uq_unstrand" # FPKM-UQ

⚠️ 定量方式怎么选：本教程的生存/预后分析统一使用 TPM；但如果要做 RNA-seq 差异表达分析（如 DESeq2），应基于原始整数 counts，不能简单把 TPM 当作 DESeq2 的输入。FPKM-UQ 是 GDC 数据体系中提供的一种标准化表达量，不同下游分析对尺度的要求不同，这里选 TPM 是本教程的分析方案，不代表所有 TCGA 预后分析都必须用 TPM。

tpm <- assay(se, "tpm_unstrand") # 行 = 基因，列 = 样本

三、单基因预后初筛：Cox + KM

3.1 开始前的分析设计

动手写代码前，先明确下面几点，能少走很多弯路：

① 分析对象：是 LUAD 还是 LUSC，还是合并的 NSCLC？

② 样本范围：是否只纳入 Primary Tumor？

③ 去重策略：如何保证”一个患者只保留一个独立样本”？

④ 生存定义：OS 的”时间”和”事件”分别怎么算？

⑤ Cox 形式：用连续表达量，还是分组变量？

⑥ KM 截断：cut-off 是否预先指定（如中位数）？

⑦ 结果定位：当前结果是”初筛/关联分析”，还是”独立外部验证”？

3.2 开始前，你需要准备什么

你只需要两样东西：

① 表达矩阵：TCGA的 SummarizedExperiment（SE）对象（含tpm_unstrand assay），或等价的”基因 × 样本” TPM 矩阵；

② 生存数据：每个样本的vital_status（Dead/Alive）、days_to_death、days_to_last_follow_up，以及样本barcode、patient。

💡 如果还没下载数据，请先按第一部分拿到SE对象。本文假设你已经有了 TCGA-LUAD_se.rds这样的文件。

📌 关键字段说明：vital_status == “Dead” 用 days_to_death 作为生存时间，否则用 ddays_to_last_follow_up（删失）。OS_event 中 1 表示死亡、0 表示删失。

3.3 核心：一个可复用脚本

下面是一个可直接套用的完整脚本。你只需要改最上面 3 行：基因名、SE 文件路径、输出目录。脚本会同时输出Cox结果、log-rank P、High vs Low 分组 Cox HR、KM 图，并保存结果表CSV。

# ================ 只需改这里 ================GENE<-"FSTL3"# ← 换成你自己的基因名SE_FILE<-"TCGA-LUAD_se.rds"# ← 换成你的 SE 文件路径OUT_DIR<-"FSTL3_validation"# ← 换成你的输出目录# =============================================library(SummarizedExperiment)library(survival)library(survminer)se<- readRDS(SE_FILE)cd<- as.data.frame(colData(se))tpm_all<- assay(se,"tpm_unstrand")# 1) 基因名 -> ENSG ID（若 SE 行名已是基因符号，可跳过）gn<- rowData(se)$gene_namenames(gn)<- rownames(se)gid<- names(gn)[gn==GENE]if(length(gid)==0)stop("未找到目标基因: ", GENE)if(length(gid)>1)warning("基因 symbol 匹配到多个 ENSG ID，取第一个: ", GENE)gid<- gid[1]# 2) 去重：1 patient = 1 sample（本数据集处理规则）# (a) 同患者、同 sample+vial+portion -> TPM 取均值（技术重复合并）# (b) 同患者多个 vial -> 保留 vial 最小 (A < B < C)，规则需在研究设计里预先声明vp<- paste(cd$patient, sapply(cd$barcode, function(bc){p<- strsplit(bc,"-")[[1]];if(length(p)>=5)paste(p[4], p[5], sep="-")elsebc}), sep="::")grp<- split(seq_len(ncol(se)), vp)tpm<- do.call(cbind, lapply(grp, function(i){if(length(i)==1)tpm_all[, i]elserowMeans(tpm_all[, i, drop=FALSE])}))clin<- do.call(rbind, lapply(grp, function(i)cd[i[1], , drop=FALSE]))# (b) 同患者多 vial -> 保留 vial 最小 (A < B < C)vl<- sapply(clin$barcode, function(bc){p<- strsplit(bc,"-")[[1]];if(length(p)>=4)substr(p[4],3,3)else"Z"})pt<- split(seq_len(nrow(clin)), clin$patient)keep<- sapply(pt, function(i)if(length(i)==1)ielsei[order(vl[i])[1]])tpm<- tpm[, keep, drop=FALSE]clin<- clin[keep, , drop=FALSE]# 3) 构建 OS 生存数据clin$OS_time<- ifelse(clin$vital_status=="Dead", as.numeric(clin$days_to_death), as.numeric(clin$days_to_last_follow_up))clin$OS_event<- ifelse(clin$vital_status=="Dead",1,0)valid<-!is.na(clin$OS_time)&clin$OS_time>0expr<- as.numeric(tpm[gid,])[valid]clin<- clin[valid, , drop=FALSE]# 4) 单因素 Cox（连续变量 log2(TPM+1)）fit_cox<- coxph(Surv(clin$OS_time, clin$OS_event)~ log2(expr +1))s<- summary(fit_cox)cox_hr<- s$conf.int[1,"exp(coef)"]cox_lo<- s$conf.int[1,"lower .95"]cox_up<- s$conf.int[1,"upper .95"]cox_p<- s$coefficients[1,"Pr(>|z|)"]cat(sprintf("Cox: HR = %.3f (%.3f-%.3f), p = %.3e\n", cox_hr, cox_lo, cox_up, cox_p))# 5) KM（中位数截断 + log-rank + High vs Low 分组 Cox HR）surv<- data.frame(time=clin$OS_time, event=clin$OS_event,expr=expr)cutoff<- median(surv$expr)surv$group<- factor(ifelse(surv$expr>cutoff,"High","Low"), levels=c("Low","High"))fit<- survfit(Surv(time, event)~ group, data=surv)lr_p<-1- pchisq(survdiff(Surv(time, event)~ group, data=surv)$chisq,df=1)grp_cox<- summary(coxph(Surv(time, event)~ group, data=surv))km_hr<- grp_cox$conf.int[1,"exp(coef)"]km_lo<- grp_cox$conf.int[1,"lower .95"]km_up<- grp_cox$conf.int[1,"upper .95"]cat(sprintf("KM: log-rank p = %.4f; High vs Low Cox HR = %.3f (%.3f-%.3f)\n", lr_p, km_hr, km_lo, km_up))# 6) 保存结果表dir.create(OUT_DIR, recursive=TRUE, showWarnings=FALSE)res<- data.frame(gene=GENE, n=nrow(clin), nevent=sum(clin$OS_event), cutoff=cutoff, cox_HR=cox_hr, cox_lower95=cox_lo, cox_upper95=cox_up, cox_p=cox_p, group_HR=km_hr, group_lower95=km_lo, group_upper95=km_up, logrank_p=lr_p)write.csv(res, file.path(OUT_DIR, paste0(GENE,"_survival_result.csv")), row.names=FALSE)# 7) 绘图p<- ggsurvplot(fit, data=surv, pval=TRUE, risk.table=TRUE, conf.int=TRUE, palette=c("#4DBBD5","#E64B35"), legend.labs=c("Low","High"), legend.title=GENE, xlab="Time (days)", ylab="Overall Survival")pdf(file.path(OUT_DIR, paste0(GENE,"_KM.pdf")), width=7, height=7)print(p)dev.off()

分步要点：

① 去重：区分”技术重复”和”独立样本”（最容易被忽略的一步）

TCGA barcode 是分层级的，例如：

TCGA-XX-XXXX-01A-01R-XXXX-XX ↑ sample+vial（p[4]=01A） ↑ portion+analyte（p[5]=01R） ↑ plate（p[6]）

生存分析通常要求一个患者只贡献一个独立样本。但”去重”不能一刀切：

• 技术重复 / 同一生物样本的重复测序：可以合并（如取均值）；

• 同一患者存在多个独立生物学样本：需要预先定义选择规则——例如优先指定样本类型、选择临床信息完整的样本，或根据研究设计决定。

🚫 避坑：本文脚本里的”同 vial 多 plate 取均值、多 vial 保留 A”是本数据集的约定处理规则，不是 TCGA 通用铁律。照搬前请确认你的数据里 A/B/C vial 确实是同一生物样本的技术分装；否则应换成你自己的选择规则，并在方法里写清楚。

② 本教程的Cox使用连续表达变量

coxph(Surv(OS_time, OS_event) ~ log2(expr + 1))

将表达量直接作为连续变量，可以减少人为 High/Low 二分造成的信息损失。但 Cox 本身并不是”只能”用连续变量——连续、二分类、多分类都可以。

⚠️ HR的单位：这里的HR对应 **log2(TPM+1)**每增加 1 个单位时的相对风险变化，而不是”原始TPM每增加 1”或严格意义的”表达翻倍”。HR > 1 表示较高表达与较高死亡风险相关，HR < 1 表示较高表达与较低死亡风险相关——这是统计学关联，不代表基因已被证明具有直接致病或保护作用。

③ KM用中位数截断，但要知道它的边界

cutoff<- median(expr)group<- ifelse(expr>cutoff,"High","Low")

中位数截断简单、透明，能避免”遍历所有截断点、选log-rank p最小的那个”这种明显过拟合。但它仍然会把连续变量硬切成两组、丢失信息；正式研究中还应结合预先定义的cut-off、连续变量模型或独立验证队列综合判断。

3.4 实战结果：以FSTL3为例

把脚本里的GENE <- “FSTL3” 跑一遍（以本次运行结果为例，经过当前样本筛选、去重和生存信息过滤后，共纳入504例患者，其中182例发生死亡事件）：

FSTL3 在 TCGA-LUAD的KM曲线

📌 结果解读：FSTL3高表达组（High，红色）的总生存低于低表达组（Low，蓝色），log-rank p = 0.0036；连续变量Cox HR = 1.29 > 1，同样指向”较高表达与较高死亡风险相关”。两者方向一致、互相印证。

⚠️ 注意：这是单因素关联分析，不能单独证明 FSTL3 是独立预后因子；正式研究还需结合临床变量做多因素Cox、检查比例风险假设，并尽可能在独立队列中验证。

📚 文献对照：本教程得到的High vs Low分组Cox HR与Meng et al. 报道的结果方向和数量级接近（文献HR = 1.55，95% CI 1.16–2.07，p = 0.003），但由于数据筛选与分析流程可能存在差异，这应理解为结果对照，而非严格意义上的完全复现。

3.5 避坑清单（汇总）

• ✅ 样本类型：GDCquery 记得 sample.type = “Primary Tumor”，避免混入其他样本类型；

• ✅ 去重：区分”技术重复”与”独立样本”，并预先声明选择规则，保证 1 患者 = 1 独立样本；

• ✅ Cox形式：本教程用连续变量 log2(TPM+1)，HR单位是”log2(TPM+1) 每增加 1”；

• ✅ KM截断：中位数截断可减少过拟合，但仍会损失连续信息；

• ✅ 比例风险假设：正式研究建议用 survival::cox.zph()检查Cox的PH假设；

• ✅ 结果定位：单因素Cox/KM是”关联/初筛”，不是”独立预后因子已被证明”；

• ✅ 换基因：仅当基因能正确匹配、且有足够表达与生存信息时，改GENE即可；否则要检查样本数、事件数与模型结果。

3.6 总结

📌 最务实的心法：第一部分拿数据，第二部分做初筛。单基因预后初筛就是”三步——把基因名一换、跑脚本、看HR和log-rank p”。同时记住：单因素 Cox/KM给出的是统计学关联，不是因果结论；连续变量Cox+合理的KM分组 + 患者级样本去重+PH假设检查，是开展单基因预后分析时值得注意的基础规范。

四、参考文献

[1]Meng X, Zhao X, Zhou B, Song W, Liang Y, Liang M, Du M, Shi J, Gao Y. FSTL3 is associated with prognosis and immune cell infiltration in lung adenocarcinoma. Journal of Cancer Research and Clinical Oncology, 2024, 150:17. DOI: 10.1007/s00432-023-05553-w