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

一、写在前面面向生信初学者的一篇实操合集分成两个独立部分第一部分教你用TCGAbiolinks从GDC下载TCGA表达数据第二部分教你拿到数据后对你自己选定的基因做单因素Cox回归 中位数截断KM生存分析以FSTL3为例把GENE - FSTL3换成自己的基因名即可。需要说明的是第二部分完成的实际上是在 TCGA 队列中的单基因预后关联分析初筛而不是严格意义上的”独立外部验证”。二、TCGA 数据下载TCGA数据获取最常用的方案是TCGAbiolinks核心四个函数查询 → 下载 → 整理 → 保存。步骤函数作用1. 查询GDCquery()按项目/数据类型/样本类型等条件生成文件清单2. 下载GDCdownload()按清单把文件下载到本地3. 整理GDCprepare()把分散文件合并成一个SummarizedExperimentSE对象4. 保存saveRDS()保存为 R 原生对象供后续分析2.1 查询GDCquerylibrary(TCGAbiolinks)library(SummarizedExperiment)query-GDCquery(projectTCGA-LUAD,# 项目 IDdata.categoryTranscriptome Profiling,# 转录组data.typeGene Expression Quantification,# 基因表达定量workflow.typeSTAR - Counts,# STAR 流程TCGA 官方推荐sample.typePrimary 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, size84.718875MB)as Sun_Sep__6_14_18_37_2026_0.tar.gz2.3 查看下载的文件结构list.files下载完成后可以用 list.files() 查看文件结构list.files(GDCdata,recursiveTRUE)# 递归列出所有文件length(list.files(GDCdata,recursiveTRUE))# 文件总数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 KM3.1 开始前的分析设计动手写代码前先明确下面几点能少走很多弯路① 分析对象是 LUAD 还是 LUSC还是合并的 NSCLC② 样本范围是否只纳入 Primary Tumor③ 去重策略如何保证”一个患者只保留一个独立样本”④ 生存定义OS 的”时间”和”事件”分别怎么算⑤ Cox 形式用连续表达量还是分组变量⑥ KM 截断cut-off 是否预先指定如中位数⑦ 结果定位当前结果是”初筛/关联分析”还是”独立外部验证”3.2 开始前你需要准备什么你只需要两样东西① 表达矩阵TCGA的 SummarizedExperimentSE对象含tpm_unstrand assay或等价的”基因 × 样本” TPM 矩阵② 生存数据每个样本的vital_statusDead/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)[gnGENE]if(length(gid)0)stop(未找到目标基因: , GENE)if(length(gid)1)warning(基因 symbol 匹配到多个 ENSG ID取第一个: , GENE)gid- gid[1]# 2) 去重1 patient 1 sample本数据集处理规则# (a) 同患者、同 samplevialportion - 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, dropFALSE])}))clin- do.call(rbind, lapply(grp, function(i)cd[i[1], , dropFALSE]))# (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)elseZ})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, dropFALSE]clin- clin[keep, , dropFALSE]# 3) 构建 OS 生存数据clin$OS_time- ifelse(clin$vital_statusDead, as.numeric(clin$days_to_death), as.numeric(clin$days_to_last_follow_up))clin$OS_event- ifelse(clin$vital_statusDead,1,0)valid-!is.na(clin$OS_time)clin$OS_time0expr- as.numeric(tpm[gid,])[valid]clin- clin[valid, , dropFALSE]# 4) 单因素 Cox连续变量 log2(TPM1)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 HRsurv- data.frame(timeclin$OS_time, eventclin$OS_event,exprexpr)cutoff- median(surv$expr)surv$group- factor(ifelse(surv$exprcutoff,High,Low), levelsc(Low,High))fit- survfit(Surv(time, event)~ group, datasurv)lr_p-1- pchisq(survdiff(Surv(time, event)~ group, datasurv)$chisq,df1)grp_cox- summary(coxph(Surv(time, event)~ group, datasurv))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, recursiveTRUE, showWarningsFALSE)res- data.frame(geneGENE, nnrow(clin), neventsum(clin$OS_event), cutoffcutoff, cox_HRcox_hr, cox_lower95cox_lo, cox_upper95cox_up, cox_pcox_p, group_HRkm_hr, group_lower95km_lo, group_upper95km_up, logrank_plr_p)write.csv(res, file.path(OUT_DIR, paste0(GENE,_survival_result.csv)), row.namesFALSE)# 7) 绘图p- ggsurvplot(fit, datasurv, pvalTRUE, risk.tableTRUE, conf.intTRUE, palettec(#4DBBD5,#E64B35), legend.labsc(Low,High), legend.titleGENE, xlabTime (days), ylabOverall Survival)pdf(file.path(OUT_DIR, paste0(GENE,_KM.pdf)), width7, height7)print(p)dev.off()分步要点① 去重区分”技术重复”和”独立样本”最容易被忽略的一步TCGA barcode 是分层级的例如TCGA-XX-XXXX-01A-01R-XXXX-XX ↑ samplevialp[4]01A ↑ portionanalytep[5]01R ↑ platep[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(TPM1)**每增加 1 个单位时的相对风险变化而不是”原始TPM每增加 1”或严格意义的”表达翻倍”。HR 1 表示较高表达与较高死亡风险相关HR 1 表示较高表达与较低死亡风险相关——这是统计学关联不代表基因已被证明具有直接致病或保护作用。③ KM用中位数截断但要知道它的边界cutoff- median(expr)group- ifelse(exprcutoff,High,Low)中位数截断简单、透明能避免”遍历所有截断点、选log-rank p最小的那个”这种明显过拟合。但它仍然会把连续变量硬切成两组、丢失信息正式研究中还应结合预先定义的cut-off、连续变量模型或独立验证队列综合判断。3.4 实战结果以FSTL3为例把脚本里的GENE - “FSTL3” 跑一遍以本次运行结果为例经过当前样本筛选、去重和生存信息过滤后共纳入504例患者其中182例发生死亡事件指标结果单因素 Cox HR (95% CI)1.29 (1.13–1.46)Cox p 值8.2 × 10⁻⁵High vs Low 分组 Cox HR (95% CI)1.54 (1.15–2.07)log-rank p0.0036中位数截断值 (TPM)34.3FSTL3 在 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.5595% CI 1.16–2.07p 0.003但由于数据筛选与分析流程可能存在差异这应理解为结果对照而非严格意义上的完全复现。3.5 避坑清单汇总✅ 样本类型GDCquery 记得 sample.type “Primary Tumor”避免混入其他样本类型✅ 去重区分”技术重复”与”独立样本”并预先声明选择规则保证 1 患者 1 独立样本✅ Cox形式本教程用连续变量 log2(TPM1)HR单位是”log2(TPM1) 每增加 1”✅ KM截断中位数截断可减少过拟合但仍会损失连续信息✅ 比例风险假设正式研究建议用 survival::cox.zph()检查Cox的PH假设✅ 结果定位单因素Cox/KM是”关联/初筛”不是”独立预后因子已被证明”✅ 换基因仅当基因能正确匹配、且有足够表达与生存信息时改GENE即可否则要检查样本数、事件数与模型结果。3.6 总结环节关键函数一句话要点数据下载GDCquery → GDCdownload → GDCprepare查、下、并拿到 SE读数据readRDSassay(se, tpm_unstrand)取 SE 的 TPM 矩阵去重两步行索引1 患者 1 独立样本生存数据vital_status days_to_*Dead 用死亡时间Alive 用随访时间单因素 Coxcoxph(Surv ~ log2(TPM1))连续变量HR 与 pKM 分析survfit survdiff ggsurvplot中位数截断 log-rank 分组 HR 最务实的心法第一部分拿数据第二部分做初筛。单基因预后初筛就是”三步——把基因名一换、跑脚本、看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