TCGA数据实战：用UCSC Xena快速搞定乳腺癌差异表达分析（附完整R代码）
在癌症研究领域，TCGA（The Cancer Genome Atlas）数据库无疑是一座金矿，它包含了33种癌症类型的基因组、转录组和表观基因组数据。但对于刚接触生物信息学的临床研究人员或学生来说，如何高效利用这些海量数据往往是个令人头疼的问题。本文将带你使用UCSC Xena平台——这个对新手极其友好的工具，快速完成乳腺癌（BRCA）数据的获取、处理和差异表达分析全流程。
不同于传统需要复杂命令行操作的TCGA数据获取方式，UCSC Xena提供了预处理的标准化数据，省去了大量数据清洗和格式转换的时间。我们将从零开始，一步步完成：
- 从Xena平台下载乳腺癌RNA-seq数据
- 使用R进行数据预处理和基因ID转换
- 利用limma包进行差异表达分析
- 生成专业级的火山图和热图可视化
所有代码都经过实战检验，可直接复制使用。即使你只有基础的R语言知识，也能在2小时内完成整个分析流程。
1. UCSC Xena平台简介与数据获取
UCSC Xena（https://xenabrowser.net/）是由加州大学圣克鲁兹分校开发的基因组数据可视化分析平台，它整合了TCGA、GTEx等多个大型项目的数据，并对原始数据进行了标准化处理。对于初学者而言，Xena有三大优势：
- 数据预处理完善：原始TCGA数据往往需要复杂的标准化和批次校正，而Xena已经完成了这些工作
- 可视化界面友好：不需要编程基础也能通过浏览器进行初步数据探索
- 下载格式统一：所有数据以相同结构存储，便于后续程序化处理
1.1 乳腺癌数据下载步骤
我们以TCGA-BRCA（乳腺癌）项目的RNA-seq数据为例：
- 访问Xena主页，选择"GDC Hub"数据集
- 在搜索框输入"BRCA"，找到"TCGA-BRCA HTSeq-FPKM"数据集
- 点击"Download"获取表达矩阵文件（通常为.tsv.gz格式）
- 同时下载基因注释文件（gencode.v22.annotation.gene.probeMap）
注意：Xena提供FPKM和counts两种表达量指标。对于差异表达分析，理论上counts更适合，但FPKM数据经过log2(fpkm+1)转换后也可使用。本文以FPKM为例。
下载完成后，你会得到两个文件：
TCGA-BRCA.htseq_fpkm.tsv.gz：表达矩阵gencode.v22.annotation.gene.probeMap：基因ID对应表
2. 数据预处理与基因ID转换
拿到原始数据后，我们需要用R进行清洗和转换。以下是完整代码及分步解释：
# 加载必要的R包 library(data.table) library(tidyverse) library(limma) # 读取表达数据（约需1-2分钟） TCGA_BRCA <- fread("TCGA-BRCA.htseq_fpkm.tsv.gz") %>% as.data.frame() # 查看数据结构 dim(TCGA_BRCA) # 行是基因，列是样本 TCGA_BRCA[1:5, 1:5] # 查看前5行5列 # 读取基因ID映射表 probeMap <- read.delim("gencode.v22.annotation.gene.probeMap") # 将Ensembl ID转换为gene symbol TCGA_BRCA <- TCGA_BRCA %>% inner_join(probeMap, by = c("Ensembl_ID" = "id")) %>% select(gene, starts_with("TCGA")) # 处理重复基因（取平均值） TCGA_BRCA <- as.data.frame( avereps(TCGA_BRCA[, -1], ID = TCGA_BRCA$gene) ) # 再次检查数据 dim(TCGA_BRCA) # 基因数量应减少 TCGA_BRCA[1:5, 1:5]2.1 关键步骤解析
- 数据导入：使用
data.table::fread()读取大型.tsv文件效率最高 - ID转换：TCGA使用Ensembl基因ID，而大多数分析需要更直观的gene symbol
- 重复基因处理：多个Ensembl ID可能对应同一gene symbol，我们取这些探针的平均值
专业提示：对于重复基因，除取平均值外，也可保留表达量最高的探针。只需将
avereps()替换为：TCGA_BRCA <- TCGA_BRCA %>% group_by(gene) %>% slice_max(rowMeans(across(starts_with("TCGA")))) %>% ungroup()
3. 样本分组与差异表达分析
TCGA样本编号包含重要信息，特别是第14-15位数字表示样本类型：
| 样本类型代码 | 含义 |
|---|---|
| 01-09 | 肿瘤组织 |
| 10-19 | 正常组织 |
3.1 构建分组信息
# 根据样本编号创建分组 TCGA_group_list <- ifelse( as.numeric(substring(colnames(TCGA_BRCA), 14, 15)) < 10, "Tumor", "Normal" ) %>% factor(levels = c("Normal", "Tumor")) # 查看各组样本数 table(TCGA_group_list)3.2 使用limma进行差异分析
虽然limma最初是为微阵列数据设计，但其voom函数使其同样适用于RNA-seq数据：
# 创建设计矩阵 design <- model.matrix(~0 + TCGA_group_list) colnames(design) <- levels(TCGA_group_list) # voom转换+线性建模 dgelist <- DGEList(counts = TCGA_BRCA, group = TCGA_group_list) %>% calcNormFactors() v <- voom(dgelist, design, plot = TRUE, normalize = "quantile") fit <- lmFit(v, design) # 设置对比组（Tumor vs Normal） contrast <- makeContrasts( contrasts = "Tumor-Normal", levels = design ) fit2 <- contrasts.fit(fit, contrast) %>% eBayes() # 提取差异表达结果 DEG <- topTable(fit2, number = Inf) %>% na.omit() DEG <- DEG[DEG$P.Value < 0.05, ] # 筛选显著差异基因 # 添加调控方向标记 DEG$regulate <- ifelse( DEG$logFC > 1 & DEG$P.Value < 0.05, "up", ifelse(DEG$logFC < -1 & DEG$P.Value < 0.05, "down", "none") ) # 查看结果 head(DEG[order(DEG$P.Value), ], 10)4. 结果可视化：火山图与热图
可视化是差异表达分析的关键环节，能直观展示整体差异模式和关键基因。
4.1 火山图绘制
火山图展示所有基因的log2倍变化与统计显著性关系：
library(ggplot2) library(ggrepel) # 标记top10上/下调基因 top_genes <- DEG %>% group_by(regulate) %>% filter(regulate %in% c("up", "down")) %>% slice_min(P.Value, n = 10) %>% pull(gene) DEG_plot <- DEG %>% mutate( label = ifelse(gene %in% top_genes, gene, NA), color = case_when( regulate == "up" ~ "#E64B35", regulate == "down" ~ "#3182BD", TRUE ~ "gray" ) ) ggplot(DEG_plot, aes(x = logFC, y = -log10(P.Value), color = color)) + geom_point(alpha = 0.6, size = 2) + geom_vline(xintercept = c(-1, 1), linetype = "dashed") + geom_hline(yintercept = -log10(0.05), linetype = "dashed") + scale_color_identity() + labs(x = "log2 Fold Change", y = "-log10 P-value") + theme_minimal(base_size = 14) + geom_text_repel( aes(label = label), max.overlaps = 20, box.padding = 0.5 )4.2 差异基因热图
热图展示关键差异基因在所有样本中的表达模式：
library(pheatmap) # 筛选显著差异基因 sig_genes <- DEG %>% filter(regulate %in% c("up", "down")) %>% pull(gene) # 提取表达矩阵（Z-score标准化） expr_matrix <- TCGA_BRCA[sig_genes, ] expr_matrix <- t(scale(t(expr_matrix))) # 样本注释信息 annotation_col <- data.frame( Group = TCGA_group_list, row.names = colnames(TCGA_BRCA) ) # 绘制热图 pheatmap( expr_matrix, annotation_col = annotation_col, show_rownames = FALSE, show_colnames = FALSE, color = colorRampPalette(c("blue", "white", "red"))(100), main = "BRCA差异表达基因热图" )5. 分析结果解读与下游分析建议
完成上述步骤后，你已获得：
- 差异表达基因列表（DEG数据框）
- 可视化图表（火山图和热图）
5.1 结果解读要点
- 火山图：每个点代表一个基因，x轴为log2倍变化（肿瘤vs正常），y轴为统计显著性。右上/左下角的点代表显著上/下调基因。
- 热图：行是差异基因，列是样本。红色表示高表达，蓝色表示低表达。观察肿瘤与正常样本是否形成明显聚类。
5.2 常见问题解决方案
问题1：差异基因数量过少
- 解决方案：调整阈值（如改用adj.P.Val代替P.Value，或放宽logFC阈值）
问题2：热图显示样本未按预期分组
- 解决方案：检查样本分组是否正确，考虑是否有批次效应需要校正
问题3：关键基因未出现在预期位置
- 解决方案：确认基因ID转换是否正确，检查原始表达量是否过低被过滤
5.3 下游分析方向
获得差异基因后，可进一步开展：
- 功能富集分析：GO、KEGG等通路分析
- 蛋白互作网络：使用STRING数据库构建关键基因网络
- 生存分析：结合临床数据评估关键基因的预后价值
- 药物靶点预测：利用CMap数据库寻找潜在治疗药物
以下是一个简单的GO富集分析代码示例：
library(clusterProfiler) library(org.Hs.eg.db) # 转换gene symbol为ENTREZ ID entrez_ids <- mapIds( org.Hs.eg.db, keys = DEG$gene[DEG$regulate == "up"], column = "ENTREZID", keytype = "SYMBOL", multiVals = "first" ) # GO富集分析 go_enrich <- enrichGO( gene = na.omit(entrez_ids), OrgDb = org.Hs.eg.db, ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.2 ) # 可视化 dotplot(go_enrich, showCategory = 20)通过本文的完整流程，即使是生物信息学新手也能在短时间内完成专业的TCGA数据分析。UCSC Xena平台大大降低了数据获取门槛，而R语言提供了灵活的分析能力。建议读者先完整复现本文流程，再根据自身研究问题调整分析策略。
