TCGA数据下载及合成基因表达矩阵


引言
在癌症研究领域,获取高质量的基因组数据对于理解疾病的本质、开发新的治疗方法至关重要。本文将详细介绍两种方法,帮助您从The Cancer Genome Atlas Program (TCGA)网站下载所需的数据,并有效管理这些宝贵资源。

以胆管癌CHOL为例,从TCGA数据库下载数据及合成基因表达矩阵。

环境准备

安装R语言及RStudio,请参考R语言安装

数据下载

方法一:直接通过TCGA官网下载数据

访问TCGA官方网站,您可以发现一个集成了大量癌症研究项目的数据库。这里不仅提供了最新的数据分析工具,如Copy Number Segment可视化工具,还定期举办网络研讨会,以指导用户更好地利用平台资源。
登录GDC Data Portal

  • 访问GDC Data Portal主页。
  • 使用您的账号登录,以便访问丰富的癌症基因组数据集。

浏览并选择项目

  • 该平台涵盖了来自86个项目的45,087例病例的高分辨率数据(截止到2025年5月7日)。
  • 您可以根据不同的原发部位(如肾、肺、乳腺等)筛选出感兴趣的研究案例。

下载数据文件
选定项目后,可以直接下载所需的临床和基因组数据文件进行进一步分析。












方法二:使用GDC客户端工具批量下载数据

为了更高效地管理和下载大规模数据集,推荐使用GDC客户端工具。

创建新文件夹datasets
在您的电脑上新建一个名为“datasets”的文件夹,用于存放即将下载的数据文件。
运行下载命令
打开终端(Mac用户请确保已移除.exe扩展名),执行如下命令:

./gdc-client download -m gdc_manifest.2025-08-04.030631.txt -d datasets

这条命令将根据指定的清单文件(例如gdc_manifest.2025-08-04.030631.txt)下载相应的数据到“datasets”目录下。请根据自己得实际情况修改保存得目录及gdc_manifest目录







无论是通过TCGA官网直接下载,还是利用GDC客户端进行批量处理,这两种方法都为研究人员提供了强大的工具来获取必要的癌症基因组数据。

合成基因表达矩阵

下载好的数据按照单个样本进行储存,我们用EXCEL随便打开一个样本的tsv观察:

TCGA新版的转录组数据把全部数据都按列打包好,放在一个文件中啦!

  • 首先,第一行告诉我们现在使用的GENECODE版本是v36(旧版是v22);
  • 其次,现在除了原本的gene_id和counts列(第一和第四列),还将gene_name、gene_type,以及tpm、fpkm、fpkm_uq等等不同标准化处理后数据一并放在了同一个表格中。

数据具体计算可参考官方:
mRNA分析流水线

表达矩阵合并

编写R脚本任意读取两个tsv文件检测gene_id和gene_name在不同样本中是否是一致的

######表达矩阵批量读入与合并######
# 先读入两个tsv检查gene_id和gene_name在不同样本中是否仍是一致的:
tsv1 <- "/0c9aba87-406f-4789-9275-e1d25bb3aea7/72fffe3e-d4fb-4862-ba01-e91289ed94ef.rna_seq.augmented_star_gene_counts.tsv"
tsv2 <- "/0e0dcae8-4cc2-4ba5-bc0d-07c9dbd0e5a2/d94cd170-9dff-4b44-a7a8-6d1c277b0f8b.rna_seq.augmented_star_gene_counts.tsv"
file1 <- paste(project_path,"/datasets/chol",tsv1,sep="")
file2 <- paste(project_path,"/datasets/chol",tsv2,sep="")
if(!file.exists(file1)){
  stop("file1不存在")
}
if(!file.exists(file2)){
  stop("file2不存在")
}
table1 <- data.table::fread(file1)
table2 <- data.table::fread(file2)
# 返回逻辑值TRUE,确认一致
gene_id_res <- identical(table1$gene_id,table2$gene_id)
if(!gene_id_res){
  stop("file1和file2的gene_id不一致")
}
# 返回逻辑值TRUE,确认一致
gene_name_res <- identical(table1$gene_name,table2$gene_name)
if(!gene_name_res){
  stop("file1和file2的gene_name不一致")
}

点击Source按钮运行脚本,如果不一致则会输出具体错误

我们选择一个任意样本,仅筛选gene_id列,用于最后表达矩阵合并完成后添加行名(因为在不同样本中gene_id的组织顺序是一致的)。

## 获取gene_id列
## 注意,data.table::fread方式读入文件并不能指定行名,所以才使用这种方法
x1 <- data.table::fread(file1,select = c("gene_id"))
head(x1)


批量读取所有tsv格式后缀文件的文件名

# 批量读取所有tsv格式后缀文件的文件名:
# *表示任意前缀,$表示固定后缀
file_list <- dir(item_path,
               pattern = "*.rna_seq.augmented_star_gene_counts.tsv$",
               recursive = T)
head(file_list)


创建自定义函数,用于批量读入所有tsv文件的unstranded列(counts)

# 创建自定义函数,用于批量读入所有tsv文件的unstranded列(counts):
exp_dt <- function(x){
  result <- data.table::fread(file.path(item_path,x),
                              select = c("unstranded"))
  return(result)
}
# 读取所有.tsv的unstranded列(将file_list转换为list,并对list中的每一个元素都应用函数exp_dt)
exp <- lapply(file_list,exp_dt)
# 将exp按列合并,并将list转化为data.table
exp <- do.call(cbind,exp)
exp[1:6,1:6]


现在表达矩阵合并好了,但还差行名(gene_id)和列名(样本名)。

# 添加行名(gene_id),所有表达矩阵的gene_id是相同的:
# data.table格式不能自定行名,因此我们先转换为数据框
exp <- as.data.frame(exp)
rownames(exp) <- x1$gene_id
exp[1:8,1:4]


现在行名有了,还差列名(TCGA样本名),列名我们需要下载对应的json文件,json文件的下载方式,可以参考上面数据下载中的第一种数据下载方式

# 列名替换为样本名
# 读取json文件(读取函数fromJSON来自R自带R包jsonlite):
json_file <- paste(json_path,"/metadata.cart.2025-08-13.json",sep="")
if(!file.exists(json_file)){
  stop("JSON文件不存在")
}
meta_dt <- jsonlite::fromJSON(json_file)
# 查看这个list中第三列中的第一个list(包含数据框):
meta_dt[[3]][[1]]


发现list中所查看的这个数据框的第一列(entity_submitter_id)即为我们所需的样本名,现在先把associated_entities这一列list组成的数据框筛选出来:

id <- meta_dt$associated_entities
id[[1]]
# 这是我们需要的样本
id[[1]][,1]
# 这是我们需要用来匹配的文件名
head(meta_dt[[4]])


批量读取所有样本ID

# 读取所有样本id,返回一个包含所有id名的向量(和前面lapply的区别,lapply返回的是list)
ids <- sapply(id,function(x){x[,1]}) 
# 生成一个包含有样本id和文件名这两列对应关系的数据框,文件名即为表达矩阵的文件名
ids2 <- data.frame(file_name = meta_dt$file_name,
                   id = ids) 
head(ids2)


观察文件名结构

ids2$file_name[1:2]
head(file_list)[1:2]

可以发现file_list中的/后面的部分也就是文件名是和file_name对应的

匹配文件名

# 用"/"做一个字符串拆分:把表达矩阵的文件名(file_list)的后半段和file_name做匹配:
library(stringr)
# 以/做拆分,simplify = T把结果返回成矩阵
file_list_dt <- str_split(file_list,"/",simplify = T)
head(file_list_dt)


检测拆分后的所有文件名,输出SUCCESS就说明全部正确

# 取矩阵第二列:
file_list_dt <- file_list_dt[,2]
# 用百分百in检测是否拆分后的文件名都完全包含在file_name中(此时顺序不同,因此不能用==)
# 全部返回TRUE,则拆分无误
res <- file_list_dt %in% ids2$file_name
if(all(res)){
  message("SUCCESS 文件名全部包含在ids的filename列中")
}else{
  stop("FAILED 文件名未全部包含在ids的filename列中")
}

调整并检查file_name顺序

# 用match函数调整file_name的顺序和file_list_dt一致(因为我们合并的表达矩阵是按照file_list_dt的顺序进行合并的)
ids3 <- ids2[match(file_list_dt,ids2$file_name),]
# 查看并检查顺序:
head(ids3$file_name)
head(file_list_dt)
res <- identical(ids3$file_name,file_list_dt)
# 返回TRUE,确认无误
if(isTRUE(res)){
  message('SUCCESS 检查通过')
}else{
  stop("FAILED 检查未通过")
}


修改列名为样本名

# 修改列名为样本名:
colnames(exp) <- ids3$id
exp[1:8,1:2]
message("表达矩阵合并完成!")



现在表达矩阵的合并我们就做好啦!

添加gene_name列

由于原始的tsv中整合了gene_name和gene_type,现在我们仅使用下载文件,同样可以快速做矩阵拆分和ID转换

## 随便选一个tsv取gene_name这一列(这里记得转换为数据框,不然后面合并会丢失exp的行名):
list <- as.data.frame(data.table::fread(file1,
                                        select = c("gene_name")))
exp2 <- cbind(list,exp)
#包含gene_name列的矩阵完成
exp2[1:10,1:3]


mRNA和lncRNA表达矩阵拆分

方法一和方法二能实现同样的结果,更推荐方法二

################# 方法一:
# 不需要区分九种不同的lncRNA了
lncRNA <- c("lncRNA")
mRNA <- c("protein_coding")

list2 <- as.data.frame(data.table::fread(file1,
                                         select = c("gene_id","gene_name","gene_type")))
lncRNA_list <- list2[list2$gene_type %in% lncRNA,]
mRNA_list <- list2[list2$gene_type %in% mRNA,]

# 现16901个,更新前14826个
length(lncRNA_list$gene_id)
# 现19962个,更新前19814个
length(mRNA_list$gene_id)

# v33版的GENECODE比v22版鉴定到的mRNA和lncRNA数量都更多了。

# 分别取表达矩阵(exp2)和mRNA/lncRNA_list的gene_id交集:
# mRNA交集:
mRNA_int <- intersect(
  rownames(exp2),
  mRNA_list$gene_id
)
# lncRNA交集:
lncRNA_int <- intersect(
  rownames(exp2),
  lncRNA_list$gene_id
)

mRNA_exp <- exp2[mRNA_int,]
lncRNA_exp <- exp2[lncRNA_int,]
###### 方法一结束

############# 方法二(直接用%in%判断并筛选即可):
# 这个方法就是我们先把表达矩阵和gene-name、gene_type三者合并,通过gene_type筛选,最后把不需要的列去掉。
list22 <- as.data.frame(data.table::fread(file1,
                                          select = c("gene_name","gene_type")))
# 把gene_type列合并进来
exp3 <- cbind(list22,exp)
exp3[1:8,1:3]

# 下面直接筛选即可:
mRNA_exp2 <- exp3[exp3$gene_type %in% c("protein_coding"),]
lncRNA_exp2 <- exp3[exp3$gene_type %in% c("lncRNA"),]

# 最后去掉用于筛选的gene_type这一列:
mRNA_exp2 <- mRNA_exp2[,-2]
lncRNA_exp2 <- lncRNA_exp2[,-2]


# 两种方法结果一致
res <- identical(mRNA_exp,mRNA_exp2)
if(!isTRUE(res)){
  stop("mRNA 方法一和方法二结果不一致")
}
identical(lncRNA_exp,lncRNA_exp2)
if(!isTRUE(res)){
  stop("lncRNA 方法一和方法二结果不一致")
}

mRNA表达矩阵:

LncRNA表达矩阵:

保存数据

# 保存数据
data_path <- paste(project_path,"/data/",item,"/",item,".Rdata",sep = "")
save(
  # 做完合并的表达矩阵
  exp,
  # 添加gene_name列的表达矩阵
  exp2,
  # mRNA的表达矩阵
  mRNA_exp,
  #lncRNA的表达矩阵
  lncRNA_exp,
  file = c(data_path)
)

至此,已经完成了TCGA数据下载及合成基因表达矩阵

脚本源代码及数据下载

脚本源代码及处理好的数据,已经开源到了GitHub上,有需要的可以直接下载获取
脚本下载

希望本文能帮助您更有效地收集和管理数据,推动科学研究向前发展

参考
TCGA更新了?!最新数据处理拿捏住

声明:初心|版权所有,违者必究|如未注明,均为原创|本网站采用BY-NC-SA协议进行授权

转载:转载请注明原文链接 - TCGA数据下载及合成基因表达矩阵


愿你勿忘初心,并从一而终