1. 为什么需要从TCGA官网获取原始count数据最近在分析TCGA-LAML数据集时我发现一个很有意思的现象从第三方获取的TPM数据竟然少了一例样本。这让我意识到依赖预处理数据可能存在潜在风险。就像做实验时要用新鲜配制的试剂一样原始数据才是分析的基础保障。count数据是RNA-seq最原始的定量结果记录了每个基因的测序reads数。而TPMTranscripts Per Million是经过长度和测序深度标准化的表达量指标。很多下游分析比如免疫浸润计算都需要TPM值但直接从TCGA下载最新count数据自行转换能确保数据完整性和可追溯性。我遇到过不少同行都有类似困扰第三方预处理数据可能存在版本滞后、样本缺失或处理参数不一致等问题。有一次合作项目中就因为使用了不同来源的TPM数据导致关键基因的表达模式出现偏差。所以现在我做分析时都会坚持从源头获取数据。2. 从TCGA官网下载count数据的完整流程2.1 准备工作首先需要安装必要的R包install.packages(c(TCGAbiolinks, SummarizedExperiment))TCGAbiolinks是TCGA官方推荐的R包就像是一个专门为TCGA数据定制的下载器。我建议创建一个专门的项目目录比如Project/ ├── Rdata/ ├── scripts/ └── ref/2.2 数据查询与下载以TCGA-LAML为例完整下载count数据的代码如下library(TCGAbiolinks) library(SummarizedExperiment) # 构建查询语句 query - GDCquery( project TCGA-LAML, data.category Transcriptome Profiling, data.type Gene Expression Quantification, workflow.type STAR - Counts ) # 下载数据约1.5GB GDCdownload(query) # 转换为R对象 dat - GDCprepare(query) laml_count - assay(dat)这里有个实用技巧可以先运行GDCquery查看返回的样本数量。我最近下载时显示151个样本与GDC官网统计一致。如果网络不稳定可以用GDCdownload(query, method client)启用断点续传。2.3 数据质量检查下载后建议立即检查# 查看前4个基因在前4个样本的表达 laml_count[1:4, 1:4] # 保存数据 save(laml_count, file Rdata/LAML_count.Rdata)重要提示新版TCGA数据使用GENCODE v36注释与老版本v22的基因ID不同。我遇到过有人混合使用不同版本导致基因匹配错误的情况。3. 获取基因有效长度的关键技术3.1 下载基因注释文件从GDC官网下载参考文件https://api.gdc.cancer.gov/data/62f23fad-0f24-43fb-8844-990d531947cf这个压缩包包含GENCODE v36注释文件。我习惯用wget下载wget -O gencode.v36.annotation.gtf.gz \ https://api.gdc.cancer.gov/data/62f23fad-0f24-43fb-8844-990d531947cf3.2 计算有效长度基因有效长度是指外显子区域的总长度。这里有个常见误区直接用CDS长度或转录本长度。实际上应该计算所有外显子的非重叠区域library(GenomicFeatures) library(parallel) # 使用多核加速节省约70%时间 cl - makeCluster(0.75*detectCores()) # 构建转录本数据库 txdb - makeTxDbFromGFF(gencode.v36.annotation.gtf.gz, formatgtf) # 按基因分组计算外显子长度 exons_gene - exonsBy(txdb, by gene) exons_gene_lens - parLapply(cl, exons_gene, function(x){ sum(width(reduce(x))) # reduce合并重叠区域 }) # 转换为数据框 geneid_efflen - data.frame( geneid names(exons_gene_lens), efflen as.numeric(exons_gene_lens) ) stopCluster(cl) save(geneid_efflen, file Rdata/gene_efflen.Rdata)关键检查点确保基因数量匹配应该是60660个基因。我有次因为内存不足导致部分基因丢失结果TPM计算全部错误。4. count到TPM的精准转换4.1 转换公式解析TPM的计算分为两步标准化基因长度标准化RPK count / (efflen/1000)测序深度标准化TPM RPK / (sum(RPK)/1e6)这就像先把每个人的收入按工作时间标准化时薪再比较相对收入水平。4.2 实际代码实现load(Rdata/gene_efflen.Rdata) load(Rdata/LAML_count.Rdata) # 必须检查基因顺序是否一致 stopifnot(identical(geneid_efflen$geneid, rownames(laml_count))) counts2TPM - function(count, efflength) { RPK - count / (efflength/1000) # 长度标准化 PMSC - sum(RPK) / 1e6 # 深度标准化因子 RPK / PMSC # TPM值 } laml_tpm - apply(laml_count, 2, counts2TPM, efflength geneid_efflen$efflen)4.3 结果验证我通常会做三个检查每列TPM总和是否为1e6允许浮点误差与已知结果对比如小洁老师的数据检查高表达基因是否合理如ACTB、GAPDHcolSums(laml_tpm)[1:5] # 应该接近1,000,0005. 常见问题与解决方案5.1 基因ID不匹配如果出现identical()返回FALSE可能是基因顺序问题# 重新排序基因 geneid_efflen - geneid_efflen[match(rownames(laml_count), geneid_efflen$geneid),]5.2 内存不足处理对于大型数据集如TCGA-GBM可以分块处理chunk_size - 10000 tpm_list - list() for(i in seq(1, nrow(laml_count), by chunk_size)){ idx - i:min(ichunk_size-1, nrow(laml_count)) tpm_list[[i]] - apply(laml_count[idx,], 2, counts2TPM, efflength geneid_efflen$efflen[idx]) } laml_tpm - do.call(rbind, tpm_list)5.3 版本兼容性问题不同GENCODE版本的基因ID系统可能不同。比如v22使用ENSG开头而v36还包含版本号ENSG00000187634.12。在整合多个数据集时要特别注意这一点。6. 完整流程自动化脚本为了方便复用我通常会封装成函数process_tcga - function(project, output_dir Rdata) { # 创建输出目录 dir.create(output_dir, showWarnings FALSE) # 下载count数据 query - GDCquery( project project, data.category Transcriptome Profiling, data.type Gene Expression Quantification, workflow.type STAR - Counts ) GDCdownload(query) dat - GDCprepare(query) count_mat - assay(dat) # 计算TPM txdb - makeTxDbFromGFF(gencode.v36.annotation.gtf.gz, formatgtf) exons_gene - exonsBy(txdb, by gene) efflen - sapply(exons_gene, function(x) sum(width(reduce(x)))) counts2TPM - function(count, efflength) { RPK - count / (efflength/1000) PMSC - sum(RPK) / 1e6 RPK / PMSC } tpm_mat - apply(count_mat, 2, counts2TPM, efflength efflen) # 保存结果 save(count_mat, file file.path(output_dir, paste0(project, _count.Rdata))) save(tpm_mat, file file.path(output_dir, paste0(project, _tpm.Rdata))) return(list(count count_mat, tpm tpm_mat)) }使用这个脚本时只需要运行result - process_tcga(TCGA-LAML)最后提醒一点TCGA数据更新时建议重新下载而不是使用本地缓存。我有次发现时隔半年后某个样本的临床注释信息被修正了。保持数据的最新状态是保证分析可靠性的基础。
