TCGA数据挖掘¶

一句话概述¶

通过GDC门户下载TCGA多组学数据，利用TCGAbiolinks进行标准化处理，开展生存分析、免疫微环境评估和突变landscape分析，是肿瘤生物信息学研究和生物标志物发现的核心数据挖掘流程。

核心知识点表格¶

知识点	关键内容	常用工具
TCGA项目概览	33种癌症、>11000样本、多组学	GDC Data Portal
数据下载	GDC API/TCGAbiolinks/gdc-client	TCGAbiolinks, GenomicDataCommons
表达分析	HTSeq counts→DESeq2/edgeR差异分析	TCGAbiolinks, DESeq2
生存分析	Kaplan-Meier/Cox回归/风险模型	survival, survminer
免疫分析	TIMER/CIBERSORT/xCell/ESTIMATE	immunedeconv, IOBR
突变分析	MAF文件解析/突变景观/驱动基因	maftools
临床关联	分子分型与临床预后的关系	survminer, forestplot
多组学整合	表达+甲基化+突变+CNV	MOFA, iClusterPlus

各步骤详解¶

第一步：理解TCGA数据资源¶

白话解释： TCGA（The Cancer Genome Atlas）是目前最大的公开癌症基因组数据库，包含33种癌症类型、超过11000个病人的多组学数据：RNA-seq表达、DNA突变、拷贝数变异、DNA甲基化、临床信息等。所有数据免费公开，是癌症研究的金矿。

技术细节： - GDC (Genomic Data Commons): TCGA数据的官方存储和分发平台 - 数据类型: RNA-seq(counts/FPKM/TPM)、WES(MAF)、SNP array(CNV)、450K/EPIC甲基化、miRNA-seq、clinical - 样本编码: TCGA-XX-XXXX-01A (01=Primary Tumor, 11=Normal Tissue) - 数据级别: Raw→Harmonized→Derived - 伦理: 公开数据无需申请，受控数据(如germline variants)需要dbGaP申请

第二步：数据下载（TCGAbiolinks）¶

代码示例：

# 安装
if (!require("BiocManager")) install.packages("BiocManager")
BiocManager::install("TCGAbiolinks")
library(TCGAbiolinks)
library(SummarizedExperiment)

# === 下载RNA-seq表达数据 ===
query_exp <- GDCquery(
    project = "TCGA-BRCA",          # 乳腺癌
    data.category = "Transcriptome Profiling",
    data.type = "Gene Expression Quantification",
    workflow.type = "STAR - Counts",
    sample.type = c("Primary Tumor", "Solid Tissue Normal"))

GDCdownload(query_exp, method = "api", files.per.chunk = 50)
data_exp <- GDCprepare(query_exp)
# 返回SummarizedExperiment对象

# 查看数据结构
dim(data_exp)  # 基因数 × 样本数
colData(data_exp)[1:5, c("patient","sample_type","vital_status","days_to_death")]

# === 下载突变数据(MAF) ===
query_maf <- GDCquery(
    project = "TCGA-BRCA",
    data.category = "Simple Nucleotide Variation",
    data.type = "Masked Somatic Mutation",
    access = "open")

GDCdownload(query_maf)
maf_data <- GDCprepare(query_maf)

# === 下载临床数据 ===
clinical <- GDCquery_clinic(project = "TCGA-BRCA", type = "clinical")
write.csv(clinical, "TCGA_BRCA_clinical.csv", row.names=FALSE)

# === 下载拷贝数数据 ===
query_cnv <- GDCquery(
    project = "TCGA-BRCA",
    data.category = "Copy Number Variation",
    data.type = "Gene Level Copy Number")

GDCdownload(query_cnv)
cnv_data <- GDCprepare(query_cnv)

第三步：差异表达分析¶

白话解释： 比较肿瘤组织和癌旁正常组织的基因表达差异，找出在癌症中上调或下调的基因。这是最基础的TCGA分析之一。

代码示例：

library(DESeq2)
library(EnhancedVolcano)

# 提取counts矩阵和样本信息
counts <- assay(data_exp, "unstranded")
coldata <- as.data.frame(colData(data_exp))

# 只保留Tumor和Normal样本
keep <- coldata$sample_type %in% c("Primary Tumor", "Solid Tissue Normal")
counts <- counts[, keep]
coldata <- coldata[keep, ]
coldata$condition <- ifelse(coldata$sample_type == "Primary Tumor", "Tumor", "Normal")

# 过滤低表达基因
keep_genes <- rowSums(counts >= 10) >= 10
counts <- counts[keep_genes, ]

# DESeq2分析
dds <- DESeqDataSetFromMatrix(
    countData = counts,
    colData = coldata,
    design = ~ condition)

dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "Tumor", "Normal"))
res_df <- as.data.frame(res[order(res$padj), ])

# 筛选DEG
sig_genes <- subset(res_df, padj < 0.01 & abs(log2FoldChange) > 2)
cat("上调:", sum(sig_genes$log2FoldChange > 0), "\n")
cat("下调:", sum(sig_genes$log2FoldChange < 0), "\n")

# 火山图
EnhancedVolcano(res_df,
    lab = rownames(res_df),
    x = 'log2FoldChange', y = 'padj',
    pCutoff = 0.01, FCcutoff = 2,
    title = 'TCGA-BRCA: Tumor vs Normal')

第四步：生存分析¶

白话解释： 生存分析回答"某个基因/分子特征是否影响病人预后"的问题。高表达某基因的病人是活得更长还是更短？通过Kaplan-Meier生存曲线和Cox回归模型来分析。

代码示例：

library(survival)
library(survminer)
library(tidyverse)

# 准备生存数据
surv_data <- clinical %>%
    select(submitter_id, vital_status, 
           days_to_death, days_to_last_follow_up) %>%
    mutate(
        time = ifelse(vital_status == "Dead", 
                      days_to_death, days_to_last_follow_up),
        status = ifelse(vital_status == "Dead", 1, 0),
        time = as.numeric(time)
    ) %>%
    filter(!is.na(time) & time > 0)

# 添加基因表达信息
# 获取标准化表达值
vsd <- vst(dds, blind=FALSE)
expr_matrix <- assay(vsd)

# 选择目标基因
gene <- "BRCA1"
gene_expr <- expr_matrix[gene, ]

# 合并
surv_data$expression <- gene_expr[match(surv_data$submitter_id, 
    substr(colnames(expr_matrix), 1, 12))]
surv_data <- surv_data[!is.na(surv_data$expression), ]

# 按中位数分高低表达组
surv_data$group <- ifelse(surv_data$expression > median(surv_data$expression),
    "High", "Low")

# Kaplan-Meier生存曲线
fit <- survfit(Surv(time, status) ~ group, data = surv_data)
p <- ggsurvplot(fit, data = surv_data,
    pval = TRUE,
    risk.table = TRUE,
    conf.int = TRUE,
    xlab = "Time (days)",
    ylab = "Overall Survival",
    title = paste(gene, "Expression and Survival"),
    palette = c("#E7B800", "#2E9FDF"))
print(p)

# Cox回归（多变量）
cox_model <- coxph(Surv(time, status) ~ expression + age + stage,
    data = surv_data)
summary(cox_model)

# 多基因风险模型（LASSO-Cox）
library(glmnet)
gene_set <- c("BRCA1","TP53","ERBB2","MKI67","ESR1","PGR")
expr_subset <- t(expr_matrix[gene_set, ])

# LASSO选择预后基因
cv_fit <- cv.glmnet(expr_subset, 
    Surv(surv_data$time, surv_data$status),
    family = "cox", alpha = 1, nfolds = 10)

# 最优lambda下的系数
coef_optimal <- coef(cv_fit, s = "lambda.min")
active_genes <- rownames(coef_optimal)[coef_optimal[,1] != 0]
cat("预后相关基因:", active_genes, "\n")

# 计算风险得分
risk_score <- predict(cv_fit, expr_subset, s="lambda.min", type="response")
surv_data$risk <- ifelse(risk_score > median(risk_score), "High_Risk", "Low_Risk")

第五步：免疫微环境评估¶

白话解释： 肿瘤免疫微环境（TIME）决定了免疫治疗的疗效。通过基因表达数据反卷积估计肿瘤中各种免疫细胞的浸润比例，评估免疫热/冷状态。

代码示例：

# 方法1: ESTIMATE算法（估计免疫/基质评分）
library(estimate)

# 准备表达矩阵(需要基因Symbol)
expr_for_estimate <- expr_matrix
# 运行ESTIMATE
filterCommonGenes(input.f="expr_matrix.gct", output.f="filtered.gct", id="GeneSymbol")
estimateScore("filtered.gct", "estimate_scores.gct")

# 方法2: 使用immunedeconv（集成多种方法）
library(immunedeconv)

# 准备TPM矩阵
tpm_matrix <- assay(data_exp, "tpm_unstrand")

# CIBERSORT (需要注册获取源文件)
# timer_result <- deconvolute(tpm_matrix, "timer", indications=rep("brca", ncol(tpm_matrix)))

# xCell
xcell_result <- deconvolute(tpm_matrix, "xcell")

# MCP-counter
mcp_result <- deconvolute(tpm_matrix, "mcp_counter")

# 方法3: ssGSEA免疫细胞评分
library(GSVA)

# 定义免疫细胞marker基因集
immune_signatures <- list(
    CD8_T = c("CD8A","CD8B","GZMA","GZMB","PRF1"),
    CD4_T = c("CD4","IL7R","CCR7","LEF1"),
    Macrophage = c("CD68","CD163","CSF1R","MSR1"),
    NK = c("NCR1","KLRD1","NKG7","GNLY"),
    B_cell = c("CD19","MS4A1","CD79A","IGHM"),
    Treg = c("FOXP3","IL2RA","CTLA4","TNFRSF18")
)

# 运行ssGSEA（GSVA >= 1.52.x 新API，旧的 gsva(..., method="ssgsea") 已弃用）
ssgsea_param <- ssgseaParam(as.matrix(tpm_matrix), immune_signatures)
immune_scores <- gsva(ssgsea_param)

# 免疫分型
library(ConsensusClusterPlus)
results <- ConsensusClusterPlus(immune_scores,
    maxK = 6, reps = 1000, pItem = 0.8,
    clusterAlg = "km", distance = "euclidean")

# 可视化免疫浸润
library(ComplexHeatmap)
Heatmap(immune_scores, 
    name = "ssGSEA score",
    show_column_names = FALSE,
    column_split = immune_clusters)

第六步：突变Landscape分析¶

白话解释： 突变景观(Mutation Landscape)全面展示癌症样本中的体细胞突变谱：哪些基因突变最频繁？什么类型的突变？与临床特征有什么关系？

代码示例：

library(maftools)

# 读取MAF文件
maf <- read.maf(maf = "TCGA-BRCA_mutations.maf.gz",
    clinicalData = clinical)

# 基本统计
getSampleSummary(maf)
getGeneSummary(maf)

# 突变景观图(Oncoplot)
oncoplot(maf, top = 20,
    clinicalFeatures = c("stage", "subtype", "vital_status"),
    sortByAnnotation = TRUE,
    removeNonMutated = FALSE)

# 突变频率柱状图
plotmafSummary(maf, rmOutlier = TRUE, addStat = 'median')

# 瀑布图（Top突变基因）
oncostrip(maf, genes = c("TP53","PIK3CA","CDH1","GATA3","MAP3K1"))

# 突变签名分析
library(MutationalPatterns)
# 或使用maftools内置
tnm <- trinucleotideMatrix(maf, ref_genome = "BSgenome.Hsapiens.UCSC.hg38")
sig <- estimateSignatures(tnm, nTry = 6)
plotSignatures(sig)

# 与COSMIC signatures比较
sig_cosmic <- compareSignatures(sig, "legacy")  # or "SBS"

# 驱动基因检测
driverGenes <- oncodrive(maf, AACol = "HGVSp_Short", minMut = 5, pvalMethod = "zscore")
plotOncodrive(driverGenes)

# 互斥/共现分析
somaticInteractions(maf, top = 25, pvalue = c(0.05, 0.01))

# 突变与生存关联
mafSurvival(maf, genes = "TP53", time = "days_to_death",
    Status = "vital_status", isTCGA = TRUE)

实战命令¶

#!/usr/bin/env Rscript
# === TCGA数据挖掘完整流程 ===

library(TCGAbiolinks)
library(DESeq2)
library(survival)
library(survminer)
library(maftools)

# 设置癌种
cancer_type <- "TCGA-LUAD"  # 肺腺癌

# 1. 下载数据
query <- GDCquery(project=cancer_type, data.category="Transcriptome Profiling",
    data.type="Gene Expression Quantification", workflow.type="STAR - Counts")
GDCdownload(query); data <- GDCprepare(query)

# 2. 差异表达
# [DESeq2分析代码]

# 3. 生存分析
clinical <- GDCquery_clinic(project=cancer_type)
# [生存分析代码]

# 4. 免疫评估
# [ESTIMATE/ssGSEA代码]

# 5. 突变分析
query_maf <- GDCquery(project=cancer_type, data.category="Simple Nucleotide Variation",
    data.type="Masked Somatic Mutation")
GDCdownload(query_maf); maf_data <- GDCprepare(query_maf)
maf <- read.maf(maf_data)
oncoplot(maf, top=20)

面试常问点¶

Q1: TCGA数据分析中，Tumor和Normal样本如何区分？ A: 通过TCGA barcode的第14-15位数字：01-09为肿瘤样本（01=Primary Tumor），10-19为正常样本（11=Solid Tissue Normal）。注意不是所有cancer type都有配对正常样本，BRCA/LUAD较多，而GBM很少。

Q2: 生存分析中，为什么不建议用最佳cutoff分组？ A: 用surv_cutpoint或X-tile找到"最佳"分组cutoff容易过拟合训练集，在独立验证中表现差。推荐使用中位数、上下四分位数等预设cutoff，或使用连续变量的Cox回归。如果必须找cutoff，应在独立验证集中确认。

Q3: TCGA RNA-seq数据用counts还是FPKM/TPM？ A: 差异表达分析必须用raw counts（DESeq2/edgeR需要原始计数进行统计检验）。样本间比较/可视化用TPM。FPKM已不推荐使用（因为不同样本间可比性差）。TCGA现在通过GDC统一使用STAR pipeline处理RNA-seq数据，GDCprepare返回的SummarizedExperiment包含unstranded、stranded_first、stranded_second、tpm_unstrand、fpkm_unstrand、fpkm_uq_unstrand多种定量，推荐差异分析用unstranded（raw counts），可视化/免疫反卷积用tpm_unstrand。

Q4: 免疫浸润分析的不同方法如何选择？ A: CIBERSORTx：金标准，可通过在线注册（https://cibersortx.stanford.edu/）免费用于学术用途，支持22种免疫细胞反卷积和单细胞参考模式。xCell：免费、覆盖64种细胞但输出是enrichment score不是比例。TIMER：简单6种免疫细胞。MCP-counter：稳健，基于标记基因。实际使用建议多种方法取交集结果更可靠。

Q5: 突变签名分析有什么生物学意义？ A: 不同突变过程（UV照射、吸烟、DNA修复缺陷、APOBEC活性等）产生特征性的突变模式（6-通道碱基替换谱）。通过反卷积可以推断肿瘤中作用的突变过程，对治疗指导有重要意义（如HRD signature预测PARP抑制剂疗效）。

易错点¶

样本配对忽略：有些分析需要配对的Tumor-Normal，需要检查patient ID匹配。
批次效应未考虑：TCGA数据来自多个中心、不同时期，存在批次效应。用DESeq2/limma时应加入batch作为协变量。
生存数据截断处理错误：存活患者的生存时间是censored（截尾数据），必须用生存分析方法（不能用t-test比较生存时间）。
多重检验未校正：同时分析数千基因与生存的关系，需要严格FDR校正。
数据版本混淆：GDC不断更新pipeline，不同时间下载的数据可能用了不同版本的处理流程。注意区分Legacy Archive（原始处理，hg19）和Harmonized Data（GDC统一pipeline，hg38）。TCGAbiolinks当前版本默认访问Harmonized数据，Legacy数据需要通过GDC Legacy Archive单独获取。
过度数据挖掘：TCGA是发现队列而非验证队列。重要发现必须在独立数据集（如GEO/ICGC/CPTAC）中验证。

补充知识¶

TCGA 33种癌症类型代码¶

缩写	全称	样本量(约)
BRCA	乳腺浸润癌	1098
LUAD	肺腺癌	585
LUSC	肺鳞癌	504
COAD	结肠腺癌	521
LIHC	肝细胞癌	377
PRAD	前列腺腺癌	498
STAD	胃腺癌	443
KIRC	肾透明细胞癌	537
GBM	胶质母细胞瘤	617
OV	卵巢浆液性癌	608

常用TCGA补充数据源¶

UCSC Xena (https://xenabrowser.net/datapages/): 预处理好的pan-cancer矩阵，下载方便，R包UCSCXenaTools可程序化下载
cBioPortal (https://www.cbioportal.org/): 在线可视化和查询
GEPIA2 (http://gepia2.cancer-pku.cn/): 在线差异表达+生存分析
TIMER2.0 (http://timer.cistrome.org/): 在线免疫浸润分析
LinkedOmics (http://www.linkedomics.org/): 在线关联分析
FireBrowse (https://firebrowse.org/): Broad提供的预计算结果

TCGA数据挖掘¶

一句话概述¶

核心知识点表格¶

各步骤详解¶

第一步：理解TCGA数据资源¶

第二步：数据下载（TCGAbiolinks）¶

第三步：差异表达分析¶

第四步：生存分析¶

第五步：免疫微环境评估¶

第六步：突变Landscape分析¶

实战命令¶

面试常问点¶

易错点¶

补充知识¶

TCGA 33种癌症类型代码¶

常用TCGA补充数据源¶

📚 相关文章推荐