做TCGA数据挖掘很多人的精力都花在差异分析本身等跑完一看结果要么差异基因数少得可怜要么多到没法解释回头排查才发现是差异分析前处理这步出了问题。这篇文章把TCGA数据处理的完整前处理流程拆开讲清楚从数据格式辨析、下载渠道选择、样本分组、表达量矩阵构建到低表达基因过滤、批次效应校正、标准化转换每个环节都说清做什么、为什么这么做、以及我踩过的坑。适合正在做TCGA数据挖掘、尤其是不太清楚数据下载之后该怎么下手的生信新手参考。1. 先搞懂TCGA数据的三种形态不然后面全是坑TCGA全称是The Cancer Genome Atlas癌症基因组图谱这是由一个国家级大规模测序项目积累下来的多组学数据库覆盖33种癌症类型、超过2万例样本包含了转录组、拷贝数变异、甲基化、突变等多层面数据。做癌症相关生信分析的人基本绕不开这个库。但很多人下载完转录组表达数据后一打开文件就懵了因为同样一个TCGA项目你能下载到好几种不同格式的表达量数据它们之间的数值含义完全不同用错了地方会让后续所有分析都失去意义。1.1 Counts、FPKM、TPM到底有什么区别TCGA的转录组表达数据在GDC官方门户和第三方平台上最常见的就是三种格式HTSeq-Counts、FPKM/FPKM-UQ、TPM。简单说HTSeq-Counts是基因比对后统计到的原始read计数它没有做任何文库大小和基因长度的校正数值大小直接受测序深度影响同一批样本之间如果测序量差异大raw counts就不适合直接拿来比较。FPKMFragments Per Kilobase Million是RNA-seq早期常用的标准化指标计算时先按基因长度归一化再按文库大小归一化。它和TPM看起来很像但计算顺序不同导致二者在数值分布上有细微差别。FPKM在跨样本比较时其实存在一点系统偏差而TPM在实现上做了修正单位是“每百万条转录本中来自该基因的转录本数”样本间可比性更好所以现在很多平台默认给TPM。那差异分析前处理到底该选哪个这取决于你后面用什么工具。如果你打算用DESeq2或edgeR跑标准差异分析官方推荐的输入是raw count矩阵HTSeq-Counts因为这两个工具内部有自己的一套归一化逻辑你给它的应该是“没被过度处理过的原始计数”。如果你只是想做基于表达量高低的分组比较、或者直接对表达值做t检验相关的差异筛选那用TPM然后log2转换是比较常见的做法。需要特别提醒绝对不要把FPKM或者TPM塞给DESeq2结果会非常难看差异基因数量动不动就成千上万明显失真。1.2 数据下载GDC官方和UCSC Xena怎么选下载TCGA表达数据主要有两条主流路径。第一条是GDC Data Portalportal.gdc.cancer.gov这是官方数据仓库能拿到最原始的HTSeq-Counts文件下载时可以选择per-sample的单独文件也可以通过GDC API或者像TCGAbiolinks这样的R包批量下载适合对数据溯源要求严格的课题。官方渠道的数据更新及时参考基因组版本和注释信息都标得很清楚写方法部分时引起来也方便。第二条是UCSC Xenaxenabrowser.net这个平台把TCGA的数据做了统一整理可以按需求下载整理好的表达矩阵提供TPM和FPKM两种格式还附带临床信息。用Xena最舒服的一点是它直接给你一个基因x样本的矩阵不需要自己逐样本合并对生信基础不太强的人友好很多。它的数据也经过了统一流程处理样本名、基因名格式都比较规范。两套方案怎么选我的建议是如果后续要用DESeq2/edgeR就走GDC拿raw count如果只做简单的表达量差异筛选Xena的TPM矩阵完全够用。另外还有一个常见渠道是cBioPortal但它的表达数据多是处理过的适合查询和可视化不太适合拿来当差异分析的原始输入。实际处理TCGA数据时经常是下载完才发现格式不对所以第一步先把数据形态和后续分析工具对应好能省下一大堆返工时间。2. 样本分组前处理的第一个分水岭拿到表达数据后的第一件事不是急着做差异分析而是把样本分组搞清楚。TCGA的样本编号不是随便起的里面藏着样本类型、组织来源、是否配对等信息如果分组出错后面所有差异比较都是白做。2.1 肿瘤和正常样本到底怎么区分TCGA样本编号是有一套规则的四段式结构比如TCGA-XX-XXXX-01A-11R-XXXX-XX其中第四段的开头两位数字代表样本类型。01开头的是原发性实体瘤11开头的是癌旁正常组织solid tissue normal05是原发性血液肿瘤06是转移性肿瘤还有02复发性肿瘤、03原发性血液肿瘤等少见类型。日常分析中最常见的就是01和11两类。这里面最容易被忽略的一点是不是所有癌症类型都同时有01和11两类样本。像GBM胶质母细胞瘤、LGG脑低级别胶质瘤这类中枢神经系统肿瘤以及一些罕见癌种正常对照样本数量极少甚至没有。这时候如果你还想着用TCGA自身样本做配对差异分析样本量就会非常尴尬。常见的替代方案是去GTEx数据库拉正常组织数据来补充对照不过这又会引入跨数据集的批次效应问题处理起来要另外花心思。在写分组文件通常叫sample_info或者phenotype表的时候有几个细节需要注意第一列名要简洁统一一般就是sample_id和group两列第二组别命名建议用明确的class标签比如Tumor和Normal而不是01和11这种数字代号后面可视化的时候图例会更清晰第三一定记得检查有没有重复的样本ID有些样本是多组织部位取样如果你下载的时候没严格过滤会出现一个病人对应多条记录的情况。2.2 表达量矩阵构建和基因ID转换如果你的数据是从GDC下载的per-sample文件每个样本是一个独立的counts文件那第一步就是把所有文件合并成一个矩阵。合并的逻辑很简单读入所有文件把里面每行基因的counts值一一对应到矩阵的列上。但实际操作有几个很容易出问题的坑。第一个坑是gene_id的格式。GDC下载的HTSeq-Counts文件里行名通常是Ensembl基因ID的带版本格式比如ENSG00000242268.11。你要用sub函数把小数点和后面的版本号去掉再去做ID注释转换。对应的基因symbol转换最常用的是用biomaRt包实时查询或者用org.Hs.eg.db包做本地注释。我自己的经验是网络不稳定时biomaRt很容易超时预先把Ensembl ID和symbol的映射关系保存成一份本地文件离线环境下也能直接转换效率高很多。第二个坑是重复基因名的处理。Ensembl ID转换到symbol之后大概率会出现多个ID对应同一个symbol的情况也就是基因名重复。这时候不能直接保留需要按累积表达量或者最大表达量合并去重否则后面差异分析会报错或者产生重复行。还有一个细节是染色体上的小RNA、假基因等非蛋白编码基因要不要过滤如果是常规mRNA表达谱差异分析我一般会保留protein_coding基因用注释文件筛掉非编码转录本能在源头减少不少噪音。3. 低表达基因过滤这一步做不做结果差非常多处理完ID转换和矩阵构建之后接下来就是过滤。别小看这一步它对后续差异分析的影响非常大但也是最容易被新手跳过的一步。3.1 为什么不能把全基因集直接拿去做差异分析TCGA的HTSeq-Counts矩阵动辄六万行但里面真正在样本中稳定表达的基因其实远没有那么多。很多基因在绝大多数样本里表达量为0或者只有几个read这些基因就是噪音。你如果不做过滤直接跑DESeq2会产生两个问题一是多重检验校正时需要比较的基因数量大幅膨胀padj会变严格一些真阳性可能被压掉二是大量全零或近零的基因会干扰后续离散度估计影响差异检验的稳健性。过滤的基本思路是设定一个保留标准比如“至少在20%的样本中counts大于等于10”这类规则。实际操作里DESeq2官方文档建议的是一个简单但有据可依的过滤方式先对所有基因计算该基因在各样本中的平均表达量或最大表达量再设一个阈值。edgeR则提供了filterByExpr函数可以综合考虑最小计数、样本量、分组信息来自动推荐过滤条件。如果你用的是TPM矩阵做后续分析通常会把阈值设成TPM 1且在多少比例样本中满足条件。3.2 过滤阈值怎么定三个场景对比我整理了几种常用的过滤标准你可以根据数据量和平时的习惯来选择。过滤方式规则示例适用场景注意事项总量阈值型行和 10 或 mean counts 5全转录组初步筛选速度快偏宽松可能留下较多低表达基因比例阈值型至少在90%样本中 counts 10关注优势表达基因的课题严格容易误删条件特异性表达基因工具推荐型filterByExpr自动判断常规差异分析准备综合库大小和分组信息个人最常用第三种是edgeR推荐型底层会综合所有样本的测序深度和库大小来设置阈值我个人用得最多因为它在保留真实信号和去掉噪音之间平衡得比较好。顺便说一句这里有一个容易混淆的概念过滤和标准化谁先谁后的问题。通常先过滤低表达基因再做标准化/归一化逻辑上更顺。原因是低表达基因的存在会影响某些标准化方法对文库大小的估计先把确定是噪音的行去掉标准化会更稳定。如果你用的是DESeq2它的median-of-ratios因素估计本身对低表达基因也敏感所以提前过滤是防患于未然。4. 批次效应处理和数据标准化批次效应是TCGA数据处理里最让人头疼的问题之一但也是差异分析前必须面对的一关。TCGA的样本不是一天之内测完的样本来源遍布多个组织中心、测序平台和批次这些技术差异如果混进你的分析里得到的结果很可能不是生物学差异而是技术噪音。4.1 批次效应怎么发现主成分分析和聚类热图批次效应是指样本在测序批次、文库制备、芯片或平台不同等因素影响下产生的系统性差异它和真实的生物学差异混杂在一起轻则让PCA图上的Tumor和Normal分不开重则直接让差异分析结果不可信。很多资料把批次效应放在差分析之后才检查但我建议在处理阶段就提前看一遍免得后续返工。最直观的方法是PCA。对log2标准化后的表达矩阵做PCA然后按样本的分组信息、测序平台如Illumina GA vs HiSeq、或者来源组织中心给点着色观察样本是否明显按非生物学因素聚类。如果发现按批次聚类的现象很明显就该考虑校正了。还有一个辅助方法是画样本相关性的热图如果同批次的样本聚成了清晰的模块而模块内既有Tumor也有Normal那基本可以判断存在批次效应。4.2 ComBat-seq和limma的removeBatchEffect怎么选处理批次效应有不少工具最常用的是sva包的ComBat系列。ComBat适合处理微阵列和RNA-seq的表达矩阵近似连续数值ComBat-seq则是针对RNA-seq raw count专门开发的它不改变数据的整数特性输出结果可以继续喂给DESeq2或edgeR做差异分析。如果你已经把数据log2转换成了连续值也可以用limma包的removeBatchEffect但它适合在标准化之后、差异分析之前对表达矩阵做残差化处理。我的组合拳实践经验是先用filterByExpr过滤掉低表达基因再用ComBat_seq对raw counts校正批次校正完再跑DESeq2比较稳健。如果你拿的是Xena的TPM矩阵TPM是连续值就log2(TPM1)之后用removeBatchEffect然后基于残差矩阵做后续分析。这里有个细节批次信息最好是样本的真实测序批次plate、seq center、tissue source site不要简单用“下载日期”替代因为下载日期跟生物学变量完全无关反而可能引入新的混淆。4.3 log2转换和标准化方法选择的经验log2转换几乎是TCGA表达数据可视化和差异筛选的标配操作但要注意两个容易被忽略的点。第一log2(x1)和log2(CPM1)的区别前者针对raw count或TPM后者是先将counts转成CPM再做log2二者数值分布不同下游算法对输入类型敏感建议从头到尾保持一致。第二log2转换只适合方差稳定的数据场景如果你要跑的是方差依赖的统计模型比如DESeq2负二项模型千万不能自己先log2再喂进去应该把raw counts原样交给DESeq2处理。至于TPM和CPM的选择我个人观点是TCGA的TPM矩阵已经考虑到基因长度影响适合做表达定量但如果你的分析目标是比较同一样本内部基因间的表达水平TPM比CPM合适。实际处理中从一个矩阵出发先明确下游分析类型再决定采用哪种数据形态不要事到临头才来回切换切换过程中数值分布的变化很可能把你的差异分析结果带偏。5. 实操从原始文件到可直接做差异分析的数据理论说了一大堆我放一套自己常用的R脚本流程出来。这个流程从GDC下载的per-sample HTSeq-Counts文件入手最终得到可以直接做差异分析的数据结构。你只需准备好两个目录一个放所有样本的counts文件每个文件两列gene_id和count一个放样本注释表格包含sample_id、group和batch等列。5.1 准备工作与读取数据library(DESeq2) library(edgeR) library(sva) library(biomaRt) library(dplyr) library(tibble) # 读入所有样本的counts文件 files - list.files(path counts_dir, pattern *.txt, full.names TRUE) sample_names - gsub(\\.txt$, , basename(files)) count_list - lapply(files, function(f) { df - read.table(f, header TRUE, row.names 1, sep \t) return(df$count) })这里有个经验如果你用GDC的下载方式文件名带长串UUID建议先重命名成样本ID避免后面矩阵列名和分组表对不上。可以用一个简单的批量重命名脚本或者直接在R里用sample_names映射总之要保持文件名和样本ID的对应关系清晰。5.2 构建表达矩阵和分组信息expr_raw - do.call(cbind, count_list) colnames(expr_raw) - sample_names # 去掉Ensembl ID的小数版本号并注释成symbol ensembl - gsub(\\..*, , rownames(expr_raw)) rownames(expr_raw) - ensembl # 用biomaRt做注释也可以预存本地映射 mart - useMart(ensembl, dataset hsapiens_gene_ensembl) annot - getBM(attributes c(ensembl_gene_id, hgnc_symbol), filters ensembl_gene_id, values ensembl, mart mart) expr_symbol - expr_raw[annot$ensembl_gene_id, ] rownames(expr_symbol) - annot$hgnc_symbol # 处理重复symbol按行求和保留 expr_symbol - expr_symbol[!is.na(rownames(expr_symbol)), ] expr_symbol - expr_symbol[rownames(expr_symbol) ! , ] expr_symbol - as.data.frame(expr_symbol) %% rownames_to_column(symbol) %% group_by(symbol) %% summarise(across(everything(), sum)) %% column_to_rownames(symbol)需要注意如果biomaRt连接不稳定建议一次性把所有Ensembl ID查完后把注释结果存成csv后面重跑时直接read.csv读取避免反复等待网络。如果你的分析不要求转symbol也可以保留Ensembl ID后续用注释文件做功能富集时再映射各有各的方便。5.3 过滤、校正、标准化三步走# 假设sample_info包含sample_id, group, batch三列 # 确保矩阵列的顺序与sample_info的sample_id完全一致 expr_raw - expr_raw[, sample_info$sample_id] # 1. 低表达基因过滤edgeR推荐方式 dge - DGEList(counts expr_raw, group sample_info$group) keep - filterByExpr(dge, group sample_info$group) dge - dge[keep, , keep.lib.sizes FALSE] # 2. 批次效应校正ComBat_seq输入raw count counts_corrected - ComBat_seq(counts dge$counts, batch sample_info$batch, group sample_info$group) # 3. 标准化转换供可视化和常规差异筛选用 # 例如转CPM后log2 expr_cpm - cpm(counts_corrected, log TRUE, prior.count 1) # 如果要喂给DESeq2做差异分析则用counts_corrected构建DESeqDataSet dds - DESeqDataSetFromMatrix(countData counts_corrected, colData sample_info, design ~ group) dds - DESeq(dds) res - results(dds, contrast c(group, Tumor, Normal))这段代码里ComBat_seq的group参数是必填的它在校正批次效应的同时会尽量保留真实的组间差异。如果你漏了这个参数ComBat_seq会在无监督模式下运行可能会把真实的生物学差异也一并“校正”掉结果就是差异分析什么都筛不出来。5.4 验证处理效果处理完之后一定要验证别急着进差异分析。常用两个检查第一重新跑一次PCA看Tumor和Normal是否按预期的分组分开了第二绘制处理前后的批次聚类热图对比确认批次效应有所缓解。如果PCA上样本仍然明显按批次聚类说明批次信息可能没找对或者批次效应与生物因素高度混杂可能需要更复杂的模型处理。PCA的可视化可以用基础R也可以ggplot2画比如提取前两个主成分按样本分组着色再用geom_text标上样本ID方便找离群点。这一步花不了五分钟但能帮你避免跑到差异分析阶段才发现数据质量有问题的大返工。6. 常见问题与排查技巧实录处理TCGA数据的过程里很多问题都是反复出现的我把一些典型场景整理成速查表遇到问题可以直接对照排查。6.1 常见报错场景速查表场景典型现象排查方向基因ID转换后全是NAbiomaRt返回大量NA检查Ensembl版本是否匹配考虑改用org.Hs.eg.db矩阵列名和分组表顺序不一致DESeq2报错样本不匹配用match()按顺序重排列彻底解决顺序问题大量基因在过滤后仍然全零过滤条件太宽松或注释比例低检查注释文件是否只覆盖了蛋白编码基因批次校正后组间差异反而变小ComBat_seq参数不当检查group参数是否正确指定不能用无监督模式TPM矩阵跑DESeq2结果高度显著但基因数异常多DESeq2不接收TPM改用TPM矩阵做线性差异筛选一个病人有多个样本同一病人在Tumor和Normal组各出现多次按病人ID去重避免伪重复混淆检验这些场景我在帮别人看代码时几乎都遇到过。尤其是矩阵列名顺序的问题看似小事跑DESeq2时一旦报错新手往往摸不着头脑其实根源就是列顺序不一致。用match函数把表达矩阵的列按sample_info的顺序重排一下问题立刻消失。6.2 避坑经验我从这些错误中学到的事最容易踩的坑是盲目照搬代码。网上的教程常常直接用Xena下载好的矩阵但你要跑的是自己从GDC下载的per-sample文件流程就不一样。我建议每一步都检查一下中间产物的行数和列数至少确认表达矩阵的基因数在过滤前后分别有多少样本数是否和分组表完全一致。数据规模对不上后面跑出什么结果都不要觉得奇怪。第二个经验是版本记录。TCGA数据本身有版本更新比如GDC上同一个TCGA项目的表达数据会随参考基因组版本更新而重新比对你下载时的release版本会直接影响Ensembl ID的注释结果。建议把下载日期、数据版本、参考基因组信息都记在一个README文件里这不仅是可重复性的要求后面写论文方法部分也会需要。第三个经验是时间成本管理。从GDC批量下载几百个per-sample文件再合并如果网速不行会比较痛苦。我曾经处理一个LUSC项目下载和整理就花了大半天后来改用TCGAbiolinks的GDCquery函数或直接从UCSC Xena拿TPM矩阵半小时内搞定。根据自己的分析目标选择合适的数据获取方式省下的时间足够你多排查好几个报错了。我做TCGA数据处理这几年最大的体会就是前处理没有想象中那么“机械”每一步都需要结合数据本身和分析目标来做决定。同样是差异分析前处理用DESeq2的人和用limma的人在过滤、标准化、批次校正的选择上可能完全不一样但核心逻辑是一致的让数据干净、可比、可解释。上面这套流程是我自己反复用过的不敢说最优但至少能帮你少走几段弯路。如果你在处理过程中遇到别的坑欢迎交流毕竟生信这条路很多经验都是踩坑踩出来的。