TCGA数据下载与整理：从GDC获取转录组数据

生信学长
2025年10月25日
TCGA数据挖掘
TCGA数据下载与整理：从GDC获取转录组数据

点击卡片关注，一起学习生信分析！



1. 背景介绍

TCGA（The Cancer Genome Atlas，癌症基因组图谱）是目前最大的癌症基因组学数据库之一，包含了33种癌症类型的多组学数据。作为生信分析的重要数据来源，掌握TCGA数据的下载和整理是进行癌症研究的基础技能。

GDC（Genomic Data Commons）数据门户是TCGA数据的官方下载平台。从GDC下载的原始数据通常需要经过整理才能用于后续分析，本文将详细介绍如何将下载的多个样本文件整理成标准的表达矩阵。

「应用场景」

  • 获取特定癌症类型的转录组数据
  • 构建基因表达矩阵用于差异分析
  • 提取样本临床信息进行生存分析

2. 分析原理

2.1 TCGA数据组织结构

从GDC下载的转录组数据通常包括：

  • 「metadata文件」 ：JSON格式，记录文件与样本的对应关系
  • 「表达数据文件」 ：每个样本一个TSV文件，包含counts、TPM、FPKM等多种表达量

2.2 数据整理流程

  1. 解析metadata文件，获取样本ID与文件名的对应关系
  2. 批量读取所有样本的表达文件
  3. 提取需要的表达量类型（如TPM）
  4. 合并成样本×基因的表达矩阵
  5. 将Ensembl ID转换为gene symbol
  6. 筛选蛋白编码基因
  7. 导出标准化的表达矩阵

3. 环境准备

3.1 R包安装 

# 安装必要的R包
install.packages("rjson")
install.packages("tidyverse")
install.packages("jsonlite")

3.2 数据准备

「第一步：从GDC门户下载数据」

  1. 访问 GDC数据门户
  2. 选择项目：TCGA-LUSC（肺鳞癌示例）
  3. 选择数据类型：Transcriptome Profiling → Gene Expression Quantification
  4. 选择分析方法：STAR - Counts
  5. 添加到购物车并下载
  6. 同时下载metadata文件（JSON格式）

「第二步：解压数据」

将下载的压缩包解压到工作目录，通常会得到一个包含多个子文件夹的目录，每个子文件夹对应一个样本。

「数据文件说明」

  • metadata.cart.2023-06-13.json ：样本元数据文件（已提供）
  • gdc_download_xxx/ ：包含所有样本数据的文件夹（需要从GDC下载）

4. 代码实现

4.1 加载R包 

# 加载必要的R包
library(rjson)
library(tidyverse)
library(jsonlite)

4.2 设置工作目录 

# 设置工作目录为数据所在文件夹
# 请根据实际情况修改路径
setwd("你的工作目录路径")

# 检查当前工作目录
getwd()

4.3 解析metadata文件

metadata文件记录了每个数据文件对应的样本信息，我们需要提取这些对应关系。 

# 读取metadata.json文件
json <- jsonlite::fromJSON("metadata.cart.2023-06-13.json")

# 查看json结构（可选）
# View(json)
# str(json)
# 提取样本ID和文件名的对应关系
sample_id <- sapply(json$associated_entities, function(x){x[,1]})

# 创建样本ID与文件名的映射表
file_sample <- data.frame(sample_id, file_name = json$file_name)

# 查看前几行
head(file_sample)

4.4 获取所有表达数据文件路径 

# 获取所有.tsv文件的路径
# 'gdc_download_xxx' 需要替换为实际的文件夹名称
count_file <- list.files('gdc_download_20230613_075345.926514/',
                         pattern = '*.tsv',
                         recursive = TRUE)

# 查看找到了多少个文件
length(count_file)

# 查看前几个文件路径
head(count_file)
# 提取文件名（不含路径）
count_file_name <- strsplit(count_file, split = '/')
count_file_name <- sapply(count_file_name, function(x){x[2]})

# 查看提取的文件名
head(count_file_name)

4.5 批量读取并合并表达数据

这是整个流程的核心步骤，我们需要逐个读取每个样本的表达文件，提取需要的列，然后合并成一个完整的表达矩阵。 

# 创建空的数据框用于存储合并后的数据
# 60660是基因数量（包括所有类型的基因）
matrix <- data.frame(matrix(nrow = 60660, ncol = 0))
# 循环读取每个样本的表达文件并合并
for (i in1:length(count_file)){
# 构建完整的文件路径
  path <- paste0('gdc_download_20230613_075345.926514//', count_file[i])

# 读取表达数据文件
  data <- read.delim(path, fill = TRUE, header = FALSE, row.names = 1)

# 第2行是列名，设置为列名
  colnames(data) <- data[2,]

# 删除前6行的注释信息
  data <- data[-c(1:6),]

# 选择需要的表达量类型（重要！）
# 列3: unstranded counts（原始counts）
# 列4: stranded_first
# 列5: stranded_second
# 列6: tpm_unstranded（TPM标准化值，推荐用于差异分析）
# 列7: fpkm_unstranded（FPKM值）
# 列8: fpkm_uq_unstranded（upper quartile标准化的FPKM）

  data <- data[6]  # 这里选择TPM值
# data <- data[3]  # 如果需要counts，取消注释这行
# data <- data[7]  # 如果需要FPKM，取消注释这行

# 将列名改为样本ID
  colnames(data) <- file_sample$sample_id[which(file_sample$file_name == count_file_name[i])]

# 合并到总表达矩阵
  matrix <- cbind(matrix, data)

# 显示进度（可选）
if(i %% 10 == 0) cat("已处理", i, "个样本\n")
}

4.6 添加基因名称和类型信息 

# 读取第一个文件获取基因名称信息
path <- paste0('gdc_download_20230613_075345.926514//', count_file[1])
data <- as.matrix(read.delim(path, fill = TRUE, header = FALSE, row.names = 1))

# 提取基因symbol
gene_name <- data[-c(1:6), 1]

# 提取基因类型
gene_type <- data[-c(1:6), 2]
# 将基因信息添加到表达矩阵
matrix0 <- cbind(gene_name, matrix)
matrix0 <- cbind(gene_type, matrix0)

# 查看矩阵维度
dim(matrix0)

4.7 数据清洗与过滤 

# 对于重复的基因symbol，保留表达量最大的那个
# 这是因为多个Ensembl ID可能对应同一个gene symbol
matrix0 <- aggregate(. ~ gene_name, data = matrix0, max)

# 检查基因数量变化
cat("合并后的基因数量:", nrow(matrix0), "\n")
# 只保留蛋白编码基因（protein_coding）
# 这一步可以根据需求选择是否执行
matrix0 <- subset(x = matrix0, gene_type == "protein_coding")

# 查看蛋白编码基因数量
cat("蛋白编码基因数量:", nrow(matrix0), "\n")

4.8 格式化并导出数据 

# 将基因名设为行名
rownames(matrix0) <- matrix0[, 1]

# 删除基因名和基因类型列
matrix0 <- matrix0[, -c(1, 2)]
# 重新整理为导出格式
matrix1 <- data.frame(ID = rownames(matrix0), matrix0)

# 将列名中的点号替换为连字符（TCGA样本ID格式）
colnames(matrix1) <- gsub('[.]', '-', colnames(matrix1))

# 查看最终矩阵的维度
dim(matrix1)

# 查看前几行和前几列
matrix1[1:5, 1:5]
# 导出为文本文件
write.table(matrix1, 'TCGA_LUSC_TPM.txt',
            sep = "\t",
            quote = FALSE,
            row.names = FALSE)

# 如果需要导出counts数据（记得修改上面的data <- data[3]）
# write.table(matrix1, 'TCGA_LUSC_count.txt',
#             sep = "\t",
#             quote = FALSE,
#             row.names = FALSE)

cat("数据导出完成！\n")

5. 结果解读

5.1 输出文件说明

运行代码后会生成 TCGA_LUSC_TPM.txt 文件，这是一个标准的基因表达矩阵：

  • 「行」 ：每一行代表一个基因（gene symbol）
  • 「列」 ：第一列是基因ID，其余列是样本ID
  • 「数值」 ：TPM标准化后的表达量

5.2 数据质量检查

「基因数量」 ：

  • 原始数据：约60,660个基因（包含所有类型）
  • 蛋白编码基因：约19,000-20,000个

「样本数量」 ：

  • 取决于下载的样本数，TCGA-LUSC通常有500+个样本
  • 包含肿瘤样本和正常样本

5.3 样本类型识别

TCGA样本ID的第14-15位标识样本类型：

  • 「01-09」 ：肿瘤样本（Tumor）
  • 「10-19」 ：正常样本（Normal）
  • 「20-29」 ：对照样本（Control）

例如： TCGA-18-3406-01A 中的 01 表示这是肿瘤样本。

5.4 表达量类型选择建议

  • 「Counts」 ：用于DESeq2等基于负二项分布的差异分析
  • 「TPM」 ：用于样本间比较，推荐用于大多数分析
  • 「FPKM」 ：较老的标准化方法，现在较少使用

6. 常见问题

6.1 文件路径错误

确保：

  • 工作目录设置正确（使用 getwd() 检查）
  • 下载的文件夹名称与代码中的一致
  • metadata.json文件在工作目录下

6.2 内存不足

处理大量样本时可能遇到内存问题，解决方法：

  • 分批处理样本
  • 增加R的内存限制： memory.limit(size = 50000)
  • 使用更高效的数据结构

6.3 基因名重复

某些Ensembl ID对应相同的gene symbol，代码中使用 aggregate() 函数保留了表达量最大的记录。也可以选择：

  • 保留平均值： aggregate(. ~ gene_name, data = matrix0, mean)
  • 保留最小值： aggregate(. ~ gene_name, data = matrix0, min)

6.4 解压文件问题

Windows系统解压可能改变文件格式，建议：


上一篇
ROC曲线分析:评估预测模型的判别能力
返回文章列表