✦
医学科研新动向
✦
TCGA
TCGA（The Cancer Genome Atlas）是一个拥有海量肿瘤基因组数据的公共数据库，为癌症研究提供了丰富的资源。掌握如何高效获取并分析这些数据对研究至关重要！
📊 本篇推文将快速了解：如何从TCGA获取肿瘤基因表达数据。
登录TCGA网站来下载数据：https://portal.gdc.cancer.gov/
/ TO /
/ TO /
01
数据下载
01. 点击 Cohort Builder
02. 以TCGA-DLBC为例：
Program → TCGA
Project → TCGA-DLBC
03. 点击Repository
04. 下载必须文件
manifest + Cart + Metadata + Clinical.csv
02
基因表达矩阵获取
将下载的文件保存在分析路径中：
# 设置工作目录，记得将其改为自己的文件路径setwd("C:\\Users\\Desktop\\data")# 指定文件名和路径metafile <- "metadata.cart.2024-10-22.json" # 下载的metadata文件名称gdcfliename <- "gdc_download_20241022_060533.962963.tar.gz" # cart文件的名称path1 <- "gdc_download_20241022_060533.962963\\" # 解压后文件的路径outfilename <- "TCGA-DLBC_FPKM.txt" # 输出表达矩阵文件的名称# 解压文件untar(gdcfliename, exdir = "gdc_download_20241022_060533.962963")# 读取metadata文件json <- jsonlite::fromJSON(metafile)sample_id <- sapply(json$associated_entities, function(x) { x[, 1] })# 创建包含样本ID和文件名的数据框file_sample <- data.frame(sample_id, file_name = json$file_name)# 列出所有包含基因计数的文件count_file <- list.files(path1, pattern = '*gene_counts.tsv', recursive = TRUE)# 提取文件名count_file_name <- strsplit(count_file, split = '/')count_file_name <- sapply(count_file_name, function(x) { x[2] })# 创建一个空矩阵，用于存储表达数据matrix <- data.frame(matrix(nrow = 60660, ncol = 0))# 循环读取每个基因计数文件，并将其合并到矩阵中for (i in 1:length(count_file_name)) {path <- paste0(path1, count_file[i]) # 构建文件路径data <- read.delim(path, fill = TRUE, header = FALSE, row.names = 1) # 读取文件colnames(data) <- data[2, ] # 设置列名为第二行内容data <- data[-c(1:6), ] # 删除前6行data <- data[7] # 选择第7列（表达数据）# 将样本ID设置为列名colnames(data) <- file_sample$sample_id[which(file_sample$file_name == count_file_name[i])]# 将数据列绑定到矩阵matrix <- cbind(matrix, data)}# 读取第一个样本的文件，用于获取基因名sample1 <- paste0(path1, count_file[1])names <- read.delim(sample1, fill = TRUE, header = FALSE, row.names = 1)colnames(names) <- names[2, ] # 设置列名为第二行内容names <- names[-c(1:6), ] # 删除前6行names <- names[, 1:2] # 选择前两列# 获取矩阵和基因名的交集same <- intersect(rownames(matrix), rownames(names))# 筛选交集部分的数据matrix <- matrix[same, ]names <- names[same, ]# 添加基因symbol到矩阵matrix$symbol <- names[, 1]matrix <- matrix[, c(ncol(matrix), 1:(ncol(matrix) - 1))] # 调整列顺序# 输出结果到指定文件write.table(matrix, file = outfilename, row.names = FALSE, quote = FALSE, sep = "\t")
03
临床数据提取
每个字段的解释：
case_submitter_id：每个样本的唯一标识符，用来标记和区分不同患者或样本。
vital_status：表示患者的生存状态，通常用“生存”（Alive）或“死亡”（Deceased）来标记。如果患者还活着，
days_to_last_follow_up会有值；如果患者去世了，days_to_death会有值。days_to_death：表示从诊断到死亡的天数。如果患者已去世，这个字段会记录他们从确诊到去世经历的时间（天数）。如果患者还活着，这个字段通常为空。
days_to_last_follow_up：从确诊到最后一次随访的天数。如果患者还活着，这个字段记录最后一次随访的时间（天数）。如果患者已去世，这个字段通常为空。
# 加载必要的库library(readr)library(tidyverse)library(dplyr)library(data.table)setwd() # clinical.tsv保存的绝对路径# 读取clinical数据clin <- fread("clinical.tsv", data.table = FALSE)# 查看数据的前几行，确保成功读取head(clin)# 提取感兴趣的列并删除重复项clin_time <- clin %>%dplyr::select(case_submitter_id,vital_status,days_to_death,days_to_last_follow_up,age_at_index,gender,ajcc_pathologic_t,ajcc_pathologic_m,ajcc_pathologic_n,ajcc_pathologic_stage) %>%dplyr::filter(!duplicated(case_submitter_id))# 查看数据的前几行以确保操作成功head(clin_time)clin_merge <- clin_time %>%dplyr::mutate(OS.time = case_when(vital_status == "Alive" ~ days_to_last_follow_up, # 若存活，使用 days_to_last_follow_upvital_status == "Dead" ~ days_to_death # 若死亡，使用 days_to_death),OS = case_when(vital_status == "Alive" ~ 0, # 存活用0表示vital_status == "Dead" ~ 1 # 死亡用1表示))# 查看数据的前几行以确保操作正确head(clin_merge)clin_merge1 <- clin_merge %>%dplyr::rename(age = age_at_index, # 将age_at_index列名修改为ageT = ajcc_pathologic_t, # 将ajcc_pathologic_t列名修改为TM = ajcc_pathologic_m, # 将ajcc_pathologic_m列名修改为MN = ajcc_pathologic_n, # 将ajcc_pathologic_n列名修改为Nstage = ajcc_pathologic_stage # 将ajcc_pathologic_stage列名修改为stage)# 查看数据的前几行以确保列名已修改head(clin_merge1)# 删除第3和第4列，并去除缺失值clin <- clin_merge1[,-(3:4)] # 删除第3和第4列clin <- na.omit(clin) # 删除包含NA的行# 将结果保存为一个制表符分隔的文件（TSV格式）write.table(clin, file = "clinical.txt", sep = "\t", row.names = FALSE, quote = FALSE)# 将结果保存为一个CSV文件write.csv(clin, file = "clinical.csv", row.names = FALSE, quote = FALSE)
-END-
文字丨本人编写，如有补充，请随时告诉我
后台回复TCGA即可获取完整代码
