R语言实战：从UCSC Xena和GDC获取TCGA/GTEx数据构建前列腺癌表达矩阵
在生物信息学研究中，TCGA和GTEx数据库是癌症基因组学研究的重要资源。本文将详细介绍如何使用R语言从UCSC Xena和GDC平台下载数据，并构建前列腺癌的TPM格式表达矩阵。这个流程不仅适用于前列腺癌研究，稍加修改后也可应用于其他癌症类型的数据分析。
1. 准备工作与环境配置
在开始数据处理前，我们需要确保R环境配置正确并安装必要的包。以下是推荐的R包列表：
# 安装必要包 install.packages(c("data.table", "jsonlite", "httr", "dplyr")) # 加载包 library(data.table) library(jsonlite) library(httr) library(dplyr)为什么选择这些包？data.table提供高效的大数据读取和处理能力；jsonlite用于解析GDC API返回的JSON数据；httr处理HTTP请求；dplyr则用于数据清洗和转换。
提示：建议使用RStudio的最新版本，并确保R版本≥4.0.0以获得最佳性能。
2. 从UCSC Xena获取GTEx数据
UCSC Xena平台整合了GTEx和TCGA等多个大型数据库的数据，提供了友好的用户界面和API接口。
2.1 数据下载
GTEx数据下载包含三个关键文件：
- 表达矩阵：log2(TPM+0.001)格式的基因表达数据
- 样本信息：包含组织来源等元数据
- 基因注释：将探针ID映射到基因名
# 设置下载URL gtex_expr_url <- "https://toil.xenahubs.net/download/gtex_RSEM_gene_tpm.gz" gtex_pheno_url <- "https://toil.xenahubs.net/download/GTEX_phenotype.gz" probe_map_url <- "https://toil.xenahubs.net/download/probeMap/gencode.v23.annotation.gene.probemap" # 下载文件 download.file(gtex_expr_url, destfile = "gtex_RSEM_gene_tpm.gz") download.file(gtex_pheno_url, destfile = "GTEX_phenotype.gz") download.file(probe_map_url, destfile = "gencode.v23.annotation.gene.probemap")2.2 数据加载与初步处理
下载完成后，我们需要将数据加载到R环境中：
# 加载表达矩阵 exp_gtex <- fread("gtex_RSEM_gene_tpm.gz", header = TRUE, sep = "\t", data.table = FALSE) rownames(exp_gtex) <- exp_gtex[,1] exp_gtex <- exp_gtex[,-1] # 加载样本信息 gtex_pheno <- fread("GTEX_phenotype.gz", header = TRUE, sep = "\t", data.table = FALSE) gtex_pheno <- gtex_pheno[,c(1,3)] # 保留样本ID和组织类型列 colnames(gtex_pheno) <- c("Barcode", "Tissue") # 筛选前列腺组织样本 prostate_samples <- gtex_pheno[gtex_pheno$Tissue == "Prostate", "Barcode"] exp_gtex <- exp_gtex[, colnames(exp_gtex) %in% prostate_samples]3. 从GDC获取TCGA数据
GDC(Genomic Data Commons)是TCGA数据的官方存储库，提供了更全面的数据但下载流程稍复杂。
3.1 通过API获取元数据
首先需要获取样本元数据，这可以通过GDC API完成：
# 设置查询参数 query <- '{ "filters": { "op": "and", "content": [ { "op": "in", "content": { "field": "cases.project.project_id", "value": ["TCGA-PRAD"] } }, { "op": "in", "content": { "field": "files.data_type", "value": ["Gene Expression Quantification"] } } ] }, "format": "JSON", "size": "1000" }' # 发送API请求 response <- POST("https://api.gdc.cancer.gov/files", body = query, encode = "json", add_headers("Content-Type" = "application/json")) # 解析响应 metadata <- fromJSON(content(response, "text"), flatten = TRUE)$data$hits3.2 下载和处理表达数据
获取元数据后，可以下载实际的表达数据文件：
# 创建下载目录 if(!dir.exists("tcga_data")) dir.create("tcga_data") # 下载函数 download_tcga_file <- function(file_id, file_name) { url <- paste0("https://api.gdc.cancer.gov/data/", file_id) dest <- paste0("tcga_data/", file_name) if(!file.exists(dest)) { GET(url, write_disk(dest, overwrite = TRUE)) } return(dest) } # 批量下载 file_paths <- mapply(download_tcga_file, metadata$file_id, metadata$file_name) # 加载并合并数据 tcga_data <- lapply(file_paths, function(x) { dt <- fread(x, data.table = FALSE) rownames(dt) <- dt[,1] return(dt[, "tpm_unstranded", drop = FALSE]) # 提取TPM列 }) exp_tcga <- do.call(cbind, tcga_data) colnames(exp_tcga) <- metadata$cases.0.submitter_id4. 数据整合与标准化
获得GTEx和TCGA数据后，需要进行一系列处理才能得到可用于分析的表达矩阵。
4.1 数据格式转换
GTEx数据需要从log2(TPM+0.001)转换回原始TPM值：
# GTEx数据转换 exp_gtex <- 2^exp_gtex - 0.001 exp_gtex[exp_gtex < 0] <- 0 # 处理可能的负值4.2 基因注释与去重
使用GENCODE注释文件对两个数据集进行基因注释：
# 加载注释文件 annot <- fread("gencode.v23.annotation.gene.probemap", header = TRUE, sep = "\t", data.table = FALSE) annot <- annot[, c("id", "gene")] # GTEx注释 common_genes <- intersect(rownames(exp_gtex), annot$id) exp_gtex <- exp_gtex[common_genes, ] gene_names <- annot[match(common_genes, annot$id), "gene"] rownames(exp_gtex) <- gene_names # TCGA注释 common_genes_tcga <- intersect(rownames(exp_tcga), annot$id) exp_tcga <- exp_tcga[common_genes_tcga, ] gene_names_tcga <- annot[match(common_genes_tcga, annot$id), "gene"] rownames(exp_tcga) <- gene_names_tcga # 去重函数 remove_duplicates <- function(expr_mat) { row_means <- rowMeans(expr_mat, na.rm = TRUE) expr_mat <- expr_mat[order(row_means, decreasing = TRUE), ] expr_mat[!duplicated(rownames(expr_mat)), ] } # 应用去重 exp_gtex <- remove_duplicates(exp_gtex) exp_tcga <- remove_duplicates(exp_tcga)4.3 数据合并
将GTEx正常组织和TCGA肿瘤组织数据合并：
# 找出共同基因 common_genes_final <- intersect(rownames(exp_gtex), rownames(exp_tcga)) # 创建合并矩阵 combined_expr <- cbind( exp_gtex[common_genes_final, ], exp_tcga[common_genes_final, ] ) # 添加样本类型信息 sample_types <- c( rep("GTEx_Normal", ncol(exp_gtex)), rep("TCGA_Tumor", ncol(exp_tcga)) ) colnames(combined_expr) <- paste0(sample_types, "_", colnames(combined_expr))5. 质量控制与数据保存
在完成表达矩阵构建后，进行必要的质量控制检查。
5.1 质量控制指标
检查数据质量的关键指标：
| 指标 | GTEx | TCGA |
|---|---|---|
| 样本数 | ncol(exp_gtex) | ncol(exp_tcga) |
| 基因数 | nrow(exp_gtex) | nrow(exp_tcga) |
| 平均表达量 | mean(rowMeans(exp_gtex)) | mean(rowMeans(exp_tcga)) |
| 零值比例 | mean(exp_gtex == 0) | mean(exp_tcga == 0) |
5.2 数据保存
将最终的表达矩阵保存为CSV和RDS格式：
# 保存为CSV write.csv(combined_expr, "prostate_cancer_expression_matrix.csv", quote = FALSE, row.names = TRUE) # 保存为RDS（保留所有信息） saveRDS(list( expression = combined_expr, sample_info = data.frame( sample_id = colnames(combined_expr), type = sample_types ), gene_info = data.frame( gene_id = rownames(combined_expr), symbol = rownames(combined_expr) ) ), file = "prostate_cancer_expression_data.rds")5.3 常见问题解决
在实际操作中可能会遇到以下问题及解决方案：
下载速度慢：
- 使用GDC Transfer Tool进行批量下载
- 考虑在云服务器上运行脚本
内存不足：
- 使用
data.table而非data.frame - 分块处理大数据文件
- 使用
基因名不一致：
- 统一使用ENSEMBL ID作为主标识符
- 考虑使用Bioconductor的注释包进行转换
样本匹配问题：
- 仔细检查元数据中的样本ID
- 使用MD5校验和验证文件完整性
