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(TPM0.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(TPM0.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 质量控制指标检查数据质量的关键指标指标GTExTCGA样本数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校验和验证文件完整性