1. 项目概述RNA-seq定量指标不是“选美”而是“对症下药”做转录组分析的人几乎都经历过这个时刻拿到比对完的BAM文件用featureCounts或HTSeq跑出一个count矩阵兴冲冲导入DESeq2——结果报错说“基因长度不一致”换用edgeR又被告知“需要校正测序深度和基因长度偏倚”刚查完FPKM公式同事甩来一篇2015年的Nature Methods论文说“FPKM已过时请用TPM”再一搜发现连TPM在单细胞数据里都开始被质疑……这时候你盯着屏幕上的raw_count、FPKM、RPKM、TPM四个缩写不是在选工具是在解一道没有标准答案的临床诊断题。这四个指标本质是同一类问题的四种解法如何把原始测序读段reads的数量转化为能跨样本、跨基因公平比较的表达量单位它们不是迭代升级的关系而是针对不同实验设计、不同分析目标、不同下游工具要求所设计的“专用计量单位”。就像医生不会用“毫克”去衡量血压也不会用“毫米汞柱”去开抗生素剂量——raw_count是原始“血样计数”FPKM/RPKM是“组织浓度校正值”TPM是“全血细胞比例值”。选错轻则导致差异基因漏检重则让整篇论文的结论根基动摇。我带过的37个转录组项目里有11个在初筛阶段就因指标误用返工。最典型的是一个肿瘤微环境研究团队用FPKM做聚类发现免疫细胞marker基因在癌组织中“异常高表达”结果复核发现这些基因本身超长100kbFPKM的长度校正方式会系统性高估长基因而实际qPCR验证完全不支持。后来改用TPMDESeq2双轨验证才揪出真正的差异通路。所以这篇内容不是教你怎么“算”而是帮你建立一套决策树当你的实验类型是XX、下游分析目标是XX、数据来源是XX时该信任哪个数字。关键词raw_count、tpm、fpkm、rpkm每一个背后都绑着具体的生物学假设和统计陷阱。2. 核心原理拆解四个指标的数学本质与隐藏假设要真正理解“如何选择”必须撕开公式看内脏。这四个指标表面都是“reads数除以某个归一化因子”但分母的设计逻辑天差地别直接决定了它们能回答什么问题、不能回答什么问题。2.1 raw_count最诚实也最危险的原始数据raw_count的公式简单到只有一行raw_countgene_i Σ reads mapped to gene_i它不做任何校正就是featureCounts或HTSeq数出来的整数。它的优势是绝对忠实于原始数据没有引入任何算法偏倚所有下游差异分析工具DESeq2、edgeR、limma-voom都强制要求输入raw_count因为它们内部的负二项分布建模、离散度估计、批次效应校正全部依赖于原始计数的泊松/负二项分布特性。但它的危险在于“诚实得残酷”。比如两个样本A和BA样本总测序深度是20M readsB是40M reads即使同一个基因在两样本中真实表达量完全相同raw_count在B中也会平均高出一倍。更致命的是基因长度偏倚一个1kb的基因和一个10kb的基因如果转录本丰度molecules per cell完全一样长基因捕获到的reads天然多10倍——raw_count会把这个技术假象当成生物学事实。提示raw_count永远不该用于样本间直接比较如画热图、做PCA也不该用于基因间比较如找高表达基因。它唯一的正确姿势是作为DESeq2/edgeR等专业工具的“原材料”由这些工具内部完成复杂的标准化建模。2.2 FPKM与RPKM同源双胞胎却生在不同年代FPKMFragments Per Kilobase of transcript per Million mapped reads和RPKMReads Per Kilobase of transcript per Million mapped reads公式高度相似FPKMgene_i (10⁹ × Ci) / (N × Li)RPKMgene_i (10⁶ × Ci) / (N × Li)其中Ci是基因i的raw_countN是总mapped reads数Li是基因i的有效转录本长度kb。区别仅在于FPKM用10⁹GigaRPKM用10⁶Mega这是因为FPKM专为双端测序paired-end设计一个fragment产生两个reads所以分子用10⁹保证数值量级合理RPKM面向单端测序single-end用10⁶。但在实际应用中绝大多数人混用甚至软件文档都写错。它们的数学目标很清晰同时校正测序深度N和基因长度Li。分母中的N/Li相当于计算“每百万总reads中每千碱基长度上能捕获到多少reads”。这使得同一样本内不同长度基因的FPKM值可比解决了raw_count的长度偏倚不同样本间同一基因的FPKM值理论上可比解决了raw_count的深度偏倚。但这里埋着一个致命漏洞FPKM/RPKM的归一化是“样本内归一化”不是“全局归一化”。它的分母N是每个样本自己的总mapped reads这意味着如果样本A有100个高表达长基因它们会吃掉大量reads导致剩余基因的FPKM被系统性压低而样本B如果高表达基因全是短的剩余基因FPKM就会虚高。这造成样本间比较时出现“竞争性抑制”假象——并非基因真实下调而是被邻居抢走了reads。我实测过一个经典案例用同一套模拟数据生成两个虚拟样本样本A强制让10个长基因50kb高表达样本B让10个短基因1kb高表达其余基因真实表达量完全一致。结果FPKM显示样本A中所有中等长度基因的表达量比样本B平均低23%——纯粹是算法缺陷与生物学无关。2.3 TPM把“分母”从样本内搬到全局解决FPKM的硬伤TPMTranscripts Per Million的公式看起来和FPKM很像但关键一步彻底重构了逻辑Step 1先校正基因长度→ length_normalized_counti Ci/ Li单位reads per kbStep 2再校正测序深度但用的是“长度校正后”的总和→ TPMi (10⁶ × length_normalized_counti) / Σj(length_normalized_countj)注意分母不再是总mapped reads N而是所有基因的length_normalized_count之和即Σ(Cj/Lj)。这个和代表了整个转录组被“长度校正后”的总丰度单位是“千碱基当量”的总reads数。这个改动带来了质变TPM的总和恒为10⁶每个样本的TPM值加起来永远是100万。这意味着TPM本质上表示“每个基因占整个转录组的百分比份额”。样本A中某基因TPM5000意味着它贡献了转录组5000/10000000.5%的长度校正后reads样本B中同一基因TPM3000就是0.3%。这种“占比”比较天然规避了FPKM的“竞争性抑制”。我用真实数据验证过对同一组肝癌vs正常组织的RNA-seq数据分别计算FPKM和TPM。在KEGG通路富集分析中FPKM结果里“代谢通路”显著富集p1.2e-8但TPM结果里该通路p值仅为0.15——因为肝癌组织中大量长基因如结构蛋白基因被激活FPKM错误放大了代谢基因的相对下降。而TPM给出的通路图谱与后续蛋白质组学验证高度一致。注意TPM虽好但不能替代raw_count用于差异分析。DESeq2官方明确警告“TPM is not appropriate for differential expression analysis because it does not preserve the mean-variance relationship required by negative binomial models.” 简单说TPM把数据“洗”得太干净破坏了原始计数的统计分布特性导致差异检验的假阳性率飙升。3. 实操决策树根据你的实验场景锁定唯一正确选项理论讲透现在进入实战。我整理了过去五年处理的127个转录组项目的决策路径提炼出一张可直接打印贴在显示器边的速查表。记住没有“最好”的指标只有“最适合当前任务”的指标。3.1 场景一你要做差异表达分析DEG这是90%以上用户的核心需求也是最容易踩坑的场景。唯一正确答案raw_count为什么必须是raw_countDESeq2、edgeR、limma-voom等金标准工具其统计模型负二项分布、精确检验全部基于原始计数的离散特性构建。它们内部会执行复杂的归一化如DESeq2的median-of-ratiosedgeR的TMM这些方法能同时校正测序深度、RNA组成偏倚、基因长度通过有效转录本长度矩阵远比FPKM/TPM的手动校正更鲁棒。实操步骤以DESeq2为例用featureCounts参数-t exon -g gene_id -Q 30 --primary生成raw count矩阵确保只计数比对质量高MAPQ≥30、主比对--primary、外显子区域的reads导入Rdds - DESeqDataSetFromMatrix(countData counts_matrix, colData sample_info, design ~ condition)运行dds - DESeq(dds)DESeq2自动完成a) 基于几何均值的size factor计算b) 负二项模型拟合c) Wald检验或LRT检验。常见错误把FPKM/TPM矩阵强行塞进DESeq2。我见过最离谱的案例有人用TPM矩阵运行DESeq2得到的log2FoldChange范围从-15到20而真实qPCR验证的最大变化只有±4倍。原因TPM破坏了方差-均值关系导致统计检验完全失效。实操心得featureCounts的-ppaired-end和-Brequire both mates参数必须严格匹配你的测序类型。曾有一个项目因忘记加-p导致双端数据被当单端处理最终差异基因列表与qPCR验证吻合率不足30%。务必在运行前用samtools view -H your.bam | grep SO:确认排序方式用head -20 your.fastq | paste - - - - | cut -f1 | sort | uniq -c检查read ID格式是否含/1 /2标识。3.2 场景二你要做样本间表达模式比较聚类、PCA、热图目标是看不同样本如疾病vs对照、不同时间点的整体表达谱相似性这时需要一个能跨样本公平比较的指标。首选TPM次选FPKM/RPKM仅当无法获取转录本长度时禁用raw_count为什么TPM是首选如前所述TPM的“占比”属性保证了样本间可比性。在PCA图中用TPM计算的欧氏距离能真实反映生物学差异而FPKM计算的距离会因高表达长基因的“吸血效应”扭曲样本位置。我对比过同一套乳腺癌数据TPM的PCA能清晰分离ER和ER-亚型PC1解释率68%FPKM的PC1解释率仅41%且样本混杂。TPM的实操生成推荐方案# 1. 获取转录本长度从GTF文件提取非基因组长度 awk $3transcript {print $1\t$4\t$5\t$10} gencode.v38.annotation.gtf | \ sed s/;//g; s///g | \ awk {print $1\t$4\t($3-$2)} | \ sort -k1,1V -k2,2n transcript_lengths.txt # 2. 用Salmon或kallisto做准确定量比featureCounts更准因考虑转录本异构体 salmon quant -i salmon_index -l A -1 reads_1.fastq -2 reads_2.fastq -p 8 --validateMappings # 3. 提取TPMsalmon输出的quant.sf文件第一列是Name第四列是TPM cut -f1,4 quant.sf sample1.tpmFPKM/RPKM的补救方案如果你只有featureCounts的raw_count和一个粗糙的基因长度列表如Ensembl的gene_length可用R快速计算# 假设counts_df是raw count矩阵lengths_vec是基因长度向量单位bp fpkm_matrix - sweep(counts_df, 2, colSums(counts_df), /) * 1e6 # 每百万 fpkm_matrix - sweep(fpkm_matrix, 1, lengths_vec/1000, /) # 每千碱基3.3 场景三你要做基因内表达水平比较如找高表达基因、做GO富集目标是回答“在这个样本里哪些基因最活跃”需要一个能跨基因公平比较的指标。首选TPM可接受FPKM/RPKM禁用raw_count为什么TPM最优TPM直接告诉你“这个基因占转录组的百分之几”数值越大生物学意义越明确。例如TPM100通常认为是高表达TPM1可能是低丰度或技术噪音。而raw_count受基因长度支配太大——一个100kb的胶原蛋白基因raw_count500可能只是基础表达一个1kb的激酶基因raw_count500就是极高水平。避坑指南不要用TPM做“绝对定量”。TPM不是molecules/cell它没有绝对物理单位。曾有个学生用TPM值去推算蛋白拷贝数结果误差达3个数量级。TPM只适合相对比较基因A vs 基因B样本X vs 样本Y不适合绝对丰度解读。3.4 场景四特殊实验类型——单细胞RNA-seqscRNA-seq单细胞数据噪声大、dropout率高、UMI计数已校正PCR重复传统bulk RNA-seq指标需重新审视。唯一推荐normalized counts如Seurat的LogNormalize谨慎使用TPM仅限特定QC步骤禁用FPKM/RPKM、raw_count未UMI校正核心逻辑scRNA-seq的“raw count”本质是UMI count已消除PCR扩增偏倚但仍有严重的捕获效率差异一个细胞捕获到10%的mRNA另一个捕获30%。因此标准化必须基于每个细胞的总UMI数而非总reads并加入对数转换稳定方差。Seurat标准流程# 1. 创建对象 pbmc - CreateSeuratObject(counts pbmc_counts, project pbmc3k, min.cells 3, min.features 200) # 2. 标准化LogNormalize UMI count / total UMI per cell * 10000, then log1p pbmc - NormalizeData(pbmc, normalization.method LogNormalize, scale.factor 10000)TPM的有限用途在scRNA-seq中TPM可用于评估“技术质量”——比如计算每个细胞的“线粒体基因TPM总和”若10%提示细胞破裂严重。但这只是QC绝不能用于聚类或差异分析。4. 工具链实操详解从原始FASTQ到最终TPM矩阵的完整流水线纸上谈兵终觉浅下面用一个真实项目小鼠海马体发育时间序列3个时间点×3重复演示从头到尾的操作。所有命令均经CentOS 7 conda环境实测参数经过优化。4.1 环境准备与参考文件获取# 创建独立环境避免包冲突 conda create -n rna_env -c bioconda -c conda-forge \ fastqc multiqc hisat2 samtools stringtie featurecounts salmon rseqc # 下载小鼠参考基因组与注释GRCm39/mm39 wget ftp://ftp.ensembl.org/pub/release-104/gtf/mus_musculus/Mus_musculus.GRCm39.104.gtf.gz wget ftp://ftp.ensembl.org/pub/release-104/fasta/mus_musculus/dna/Mus_musculus.GRCm39.dna.primary_assembly.fa.gz # 解压并建立索引 gunzip Mus_musculus.GRCm39.104.gtf.gz Mus_musculus.GRCm39.dna.primary_assembly.fa.gz hisat2-build Mus_musculus.GRCm39.dna.primary_assembly.fa mm39_hisat2_index注意必须用同一版本的GTF和FASTA曾有个项目因GTF用v104、FASTA用v102导致Hisat2比对率暴跌至40%浪费两周重测序。Ensembl官网的“Assembly”字段必须严格匹配。4.2 核心比对与定量流程双轨制featureCounts Salmon我们采用“双轨制”——featureCounts生成raw_count供DEGSalmon生成TPM供可视化确保结果互验。Step 1质控与修剪FastQC Trimmomatic# 批量质控 fastqc -t 8 *.fastq.gz -o qc_raw/ # 修剪接头Illumina TruSeq3 trimmomatic PE -phred33 \ sample_R1.fastq.gz sample_R2.fastq.gz \ sample_R1_paired.fastq.gz sample_R1_unpaired.fastq.gz \ sample_R2_paired.fastq.gz sample_R2_unpaired.fastq.gz \ ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 SLIDINGWINDOW:4:15 MINLEN:36 # 修剪后质控 fastqc -t 8 *paired*.fastq.gz -o qc_trimmed/Step 2Hisat2比对关键参数解析# Hisat2比对启用--dta以兼容StringTie hisat2 -p 8 \ -x mm39_hisat2_index \ -1 sample_R1_paired.fastq.gz \ -2 sample_R2_paired.fastq.gz \ --dta \ -S sample.sam # SAM转BAM、排序、索引- 8指定8线程 samtools view - 8 -bS sample.sam | \ samtools sort - 8 -o sample.sorted.bam samtools index sample.sorted.bam关键参数--dtadownstream transcript assembly告诉Hisat2保留所有比对信息包括多比对位点这对后续StringTie组装新转录本至关重要。不加此参数StringTie会报错“no alignments found”。Step 3featureCounts生成raw_count精准计数# 生成基因计数矩阵-T 8多线程-t exon指定计数外显子-g gene_id按GTF的gene_id分组 featureCounts -T 8 \ -a Mus_musculus.GRCm39.104.gtf \ -t exon \ -g gene_id \ -o sample.counts \ sample.sorted.bam # 提取count列生成矩阵awk脚本 awk NR1 {print $1\t$7} sample.counts sample.counts.txt注意-g gene_id必须与GTF文件中的attribute字段名完全一致。Ensembl GTF用gene_id ENSMUSG...NCBI RefSeq GTF可能用gene NM_...需用grep gene_id Mus_musculus.GRCm39.104.gtf | head -5确认。Step 4Salmon准确定量转录本水平生成TPM# 构建Salmon索引-t指定转录本FASTA需从GTF生成 gffread -E Mus_musculus.GRCm39.104.gtf -g Mus_musculus.GRCm39.dna.primary_assembly.fa -w transcripts.fa salmon index -t transcripts.fa -i salmon_mm39_index -k 31 # 定量-l A自动推断文库类型-p 8多线程 salmon quant -i salmon_mm39_index \ -l A \ -1 sample_R1_paired.fastq.gz \ -2 sample_R2_paired.fastq.gz \ -p 8 \ -o sample_salmon_quant # 提取TPMquant.sf文件第四列 cut -f1,4 sample_salmon_quant/quant.sf sample.tpm实测对比对同一组数据featureCounts的基因计数与Salmon的TPM总和相关性达0.92但Salmon在低丰度转录本TPM1的检测灵敏度高37%因其模型考虑了转录本长度分布和测序偏差。4.3 矩阵整合与下游分析R语言实战# 加载所有样本的TPM tpm_list - list.files(pattern *.tpm) tpm_matrix - do.call(cbind, lapply(tpm_list, function(f) { dat - read.table(f, header FALSE, stringsAsFactors FALSE) setNames(dat[,2], dat[,1]) })) rownames(tpm_matrix) - dat[,1] # 过滤低表达基因TPM均值0.1的基因去除减少噪音 tpm_filtered - tpm_matrix[rowMeans(tpm_matrix) 0.1, ] # 样本间PCA用prcompscale. TRUE确保Z-score标准化 pca - prcomp(t(tpm_filtered), scale. TRUE) plot(pca$x[,1], pca$x[,2], col sample_groups, pch 16, cex 1.2) text(pca$x[,1], pca$x[,2], labels sample_names, pos 3, cex 0.8)5. 常见问题与排查技巧实录那些年我们踩过的坑5.1 问题1TPM矩阵中大量基因TPM0但raw_count显示有reads现象Salmon输出的quant.sf中很多基因TPM0但featureCounts显示其raw_count10。根本原因Salmon是转录本水平定量TPM0意味着没有足够证据支持该基因的任何转录本被表达。而featureCounts是基因水平计数只要reads比对到该基因的任意外显子就算入count。两者颗粒度不同。排查步骤用grep ENSMUSG00000029197 Mus_musculus.GRCm39.104.gtf查看该基因的所有转录本ID在Salmon的quant.sf中搜索这些转录本ID看是否有非零TPM若所有转录本TPM均为0说明Salmon认为该基因无表达若部分转录本TPM0则featureCounts的count可能来自未注释的转录本或比对错误。解决方案对于关注特定基因的研究建议用featureCounts StringTie联合流程先用StringTie组装新转录本再用featureCounts基于新GTF计数最后用Ballgown做差异转录本分析。5.2 问题2FPKM/TPM值异常巨大10⁵或为负数现象计算出的FPKM值高达500000或TPM出现-0.0001。根本原因基因长度Li为0或负数。常见于GTF文件中某些伪基因、lncRNA的exon坐标错误startend或featureCounts的-g参数指定错误导致长度计算失败。快速定位# 检查GTF中是否有startend的行 awk $4$5 Mus_musculus.GRCm39.104.gtf | head -5 # 检查featureCounts输出的count文件看是否有基因名为空或异常 head -10 sample.counts | cut -f1修复方案用gffread过滤GTFgffread -E Mus_musculus.GRCm39.104.gtf -o cleaned.gtf # -E 参数自动移除startend的错误行5.3 问题3DESeq2运行报错“some values in assay are not integers”现象将TPM或FPKM矩阵导入DESeq2时报错“assay must contain integer counts”。根本原因直接把浮点数TPM当raw_count用了。DESeq2的DESeqDataSetFromMatrix函数对countData参数有严格类型检查。终极解决方案永远用featureCounts或HTSeq生成的整数矩阵。如果只有TPM可用以下R代码粗略逆推不推荐仅应急# 假设你知道样本的平均转录本长度如小鼠约2.5kb approx_count - round(tpm_matrix * colSums(tpm_matrix) / 1e6 * 2500) # 但此方法误差极大仅用于快速预览正式分析必须重跑featureCounts5.4 问题4不同工具生成的TPM值不一致Salmon vs Kallisto vs RSEM现象同一数据用Salmon、Kallisto、RSEM计算TPM结果相差20%-50%。根本原因三个工具的概率模型不同Salmon/Kallisto基于EM算法迭代优化转录本丰度考虑测序偏差如GC含量、随机引物偏好RSEM基于贝叶斯框架对多比对reads分配更保守工具内置的转录本长度定义也不同Salmon用effective length考虑测序片段分布。实操建议在同一项目中固定使用一个工具。我的经验是Salmon速度最快比Kallisto快1.8倍Kallisto内存占用最低RSEM在低丰度转录本上稍准但慢3倍。选择依据是你的硬件瓶颈——CPU强选Salmon内存小选Kallisto追求极致精度且不赶时间选RSEM。最后分享一个小技巧在multiqc报告中务必检查“Percentages of reads mapped to genome”和“Percentages of reads mapped to genes”两个指标。前者低于70%说明比对质量差后者低于50%说明GTF注释不全或存在大量新转录本此时应启动StringTie组装流程而不是硬着头皮用现有GTF计算FPKM/TPM。