5步自动化流程：用R语言高效处理TCGA食管癌RNA-seq数据
每次拿到TCGA的原始数据，你是不是也经历过这样的痛苦？解压后几十个样本文件夹散落各处，metadata文件像天书一样难懂，基因ID转换总出问题，手动合并数据耗时又容易出错。作为生物信息学分析的基础环节，数据清洗的质量直接影响后续差异分析、生存分析的可靠性。今天我们就用R语言打造一个全自动流水线，从原始文件到清洗后的表达矩阵一气呵成。
1. 环境准备与数据解压
工欲善其事，必先利其器。在开始前需要确保以下R包已安装：
# 安装必要包（若未安装） install.packages(c("tidyverse", "rjson", "data.table", "biomaRt"))TCGA数据通常以Cart形式下载，包含两个关键文件：
gdc_manifest.txt：记录所有样本文件的哈希值和路径metadata.json：包含样本与文件对应关系的元数据
解压后目录结构通常如下：
TCGA-ESCA/ ├── 01f51a3d-7c04-4b2a-9e5f-1f3e8a2b4c5c/ │ ├── 01f51a3d-7c04-4b2a-9e5f-1f3e8a2b4c5c.rna_seq.augmented_star_gene_counts.tsv ├── 0a2b4c5c-7c04-4b2a-9e5f-1f51a3d7c04b/ │ ├── 0a2b4c5c-7c04-4b2a-9e5f-1f51a3d7c04b.rna_seq.augmented_star_gene_counts.tsv ...提示：建议使用
list.files(recursive = TRUE)快速查看解压后的文件结构，确认所有样本文件均包含.tsv后缀的表达量数据。
2. 元数据解析与样本匹配
元数据是连接文件名与样本ID的桥梁。以下代码自动提取样本信息：
library(tidyverse) library(rjson) # 读取元数据 meta <- fromJSON(file = "metadata.json") # 提取样本信息 samples <- map_dfr(meta, ~{ tibble( file_id = .x$file_id, sample_id = .x$associated_entities[[1]]$entity_submitter_id, file_name = .x$file_name ) }) # 显示前5个样本 head(samples, 5)输出示例：
| file_id | sample_id | file_name |
|---|---|---|
| 01f51a3d-7c04-4b2a-9e5f-1f3e8a2b4c5c | TCGA-AA-1234-01 | 01f51a3d...augmented_star_gene_counts.tsv |
| 0a2b4c5c-7c04-4b2a-9e5f-1f51a3d7c04b | TCGA-BB-5678-11 | 0a2b4c5c...augmented_star_gene_counts.tsv |
关键点解析：
entity_submitter_id中的"-01"表示原发肿瘤，"-11"表示正常组织- 文件名中的UUID需要与
gdc_manifest.txt中的记录对应
3. 批量读取与表达矩阵构建
使用data.table包高效读取多个文件：
library(data.table) # 获取所有数据文件路径 file_paths <- list.files(pattern = "*.tsv", recursive = TRUE, full.names = TRUE) # 自定义读取函数 read_tsv_fast <- function(path) { dt <- fread(path, select = c("gene_id", "unstranded")) setnames(dt, "unstranded", basename(path)) return(dt) } # 并行读取（加快速度） expr_list <- lapply(file_paths, read_tsv_fast) # 合并所有样本 expr_matrix <- reduce(expr_list, full_join, by = "gene_id")注意：内存不足时可分批次读取，每次处理10-20个样本后再合并
4. 基因注释与过滤
ENSEMBL ID需要转换为更易读的Gene Symbol：
library(biomaRt) # 连接ENSEMBL数据库 mart <- useMart("ensembl", dataset = "hsapiens_gene_ensembl") # 获取基因注释 annot <- getBM( attributes = c("ensembl_gene_id", "hgnc_symbol", "gene_biotype"), filters = "ensembl_gene_id", values = expr_matrix$gene_id, mart = mart ) # 合并注释信息 expr_annot <- expr_matrix %>% left_join(annot, by = c("gene_id" = "ensembl_gene_id")) %>% filter(gene_biotype == "protein_coding") %>% select(-gene_id, -gene_biotype) %>% group_by(hgnc_symbol) %>% summarise(across(where(is.numeric), mean)) %>% column_to_rownames("hgnc_symbol")常见问题处理：
- 重复基因名：取表达量平均值（也可选择最大值）
- 非编码RNA：通过
gene_biotype过滤掉"pseudogene"等类型
5. 质量控制与数据导出
最后的质量控制步骤：
# 过滤低表达基因 keep <- rowSums(expr_annot > 0) >= 0.5*ncol(expr_annot) expr_final <- expr_annot[keep, ] # 样本聚类检查离群值 hc <- hclust(dist(t(log2(expr_final+1)))) plot(hc, main = "Sample Clustering") # 保存结果 write.csv(expr_final, "TCGA_ESCA_Counts_Matrix.csv")完整流程不到50行代码，但实现了：
- 自动样本匹配
- 并行文件读取
- 智能基因注释
- 一键式质控
这个流程同样适用于其他TCGA癌症类型，只需替换元数据文件和样本路径即可。我曾用这套方法处理过TCGA-LUAD的500+样本，原本需要一整天的手工操作现在只需15分钟就能完成。
