企业官网建设流程全解析

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基因在癌组织中“异常高表达”，结果复核发现，这些基因本身超长（>100kb），FPKM的长度校正方式会系统性高估长基因，而实际qPCR验证完全不支持。后来改用TPM+DESeq2双轨验证，才揪出真正的差异通路。所以这篇内容不是教你怎么“算”，而是帮你建立一套决策树：当你的实验类型是XX、下游分析目标是XX、数据来源是XX时，该信任哪个数字。关键词raw_count、tpm、fpkm、rpkm，每一个背后都绑着具体的生物学假设和统计陷阱。

2. 核心原理拆解：四个指标的数学本质与隐藏假设

要真正理解“如何选择”，必须撕开公式看内脏。这四个指标表面都是“reads数除以某个归一化因子”，但分母的设计逻辑天差地别，直接决定了它们能回答什么问题、不能回答什么问题。

2.1 raw_count：最诚实也最危险的原始数据

raw_count的公式简单到只有一行：
raw_count_{gene_i}= Σ reads mapped to gene_i

它不做任何校正，就是featureCounts或HTSeq数出来的整数。它的优势是绝对忠实于原始数据：没有引入任何算法偏倚，所有下游差异分析工具（DESeq2、edgeR、limma-voom）都强制要求输入raw_count，因为它们内部的负二项分布建模、离散度估计、批次效应校正，全部依赖于原始计数的泊松/负二项分布特性。

但它的危险在于“诚实得残酷”。比如两个样本A和B，A样本总测序深度是20M reads，B是40M reads，即使同一个基因在两样本中真实表达量完全相同，raw_count在B中也会平均高出一倍。更致命的是基因长度偏倚：一个1kb的基因和一个10kb的基因，如果转录本丰度（molecules per cell）完全一样，长基因捕获到的reads天然多10倍——raw_count会把这个技术假象当成生物学事实。

提示：raw_count永远不该用于样本间直接比较（如画热图、做PCA），也不该用于基因间比较（如找高表达基因）。它唯一的正确姿势是作为DESeq2/edgeR等专业工具的“原材料”，由这些工具内部完成复杂的标准化建模。

2.2 FPKM与RPKM：同源双胞胎，却生在不同年代

FPKM（Fragments Per Kilobase of transcript per Million mapped reads）和RPKM（Reads Per Kilobase of transcript per Million mapped reads）公式高度相似：

FPKM_{gene_i}= (10⁹ × C_i) / (N × L_i)
RPKM_{gene_i}= (10⁶ × C_i) / (N × L_i)

其中C_i是基因i的raw_count，N是总mapped reads数，L_i是基因i的有效转录本长度（kb）。区别仅在于FPKM用10⁹（Giga），RPKM用10⁶（Mega），这是因为FPKM专为双端测序（paired-end）设计，一个fragment产生两个reads，所以分子用10⁹保证数值量级合理；RPKM面向单端测序（single-end），用10⁶。但在实际应用中，绝大多数人混用，甚至软件文档都写错。

它们的数学目标很清晰：同时校正测序深度（N）和基因长度（L_i）。分母中的N/L_i相当于计算“每百万总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的硬伤

TPM（Transcripts Per Million）的公式看起来和FPKM很像，但关键一步彻底重构了逻辑：

Step 1：先校正基因长度
→ length_normalized_count_i= C_i/ L_i（单位：reads per kb）

Step 2：再校正测序深度（但用的是“长度校正后”的总和）
→ TPM_i= (10⁶ × length_normalized_count_i) / Σ_j(length_normalized_count_j)

注意分母：不再是总mapped reads N，而是所有基因的length_normalized_count之和，即Σ(C_j/L_j)。这个和代表了整个转录组被“长度校正后”的总丰度，单位是“千碱基当量”的总reads数。

这个改动带来了质变：

• TPM的总和恒为10⁶：每个样本的TPM值加起来永远是100万。这意味着TPM本质上表示“每个基因占整个转录组的百分比份额”。样本A中某基因TPM=5000，意味着它贡献了转录组5000/1000000=0.5%的长度校正后reads；样本B中同一基因TPM=3000，就是0.3%。这种“占比”比较，天然规避了FPKM的“竞争性抑制”。

我用真实数据验证过：对同一组肝癌vs正常组织的RNA-seq数据，分别计算FPKM和TPM。在KEGG通路富集分析中，FPKM结果里“代谢通路”显著富集（p=1.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_count？DESeq2、edgeR、limma-voom等金标准工具，其统计模型（负二项分布、精确检验）全部基于原始计数的离散特性构建。它们内部会执行复杂的归一化（如DESeq2的median-of-ratios，edgeR的TMM），这些方法能同时校正测序深度、RNA组成偏倚、基因长度（通过有效转录本长度矩阵），远比FPKM/TPM的手动校正更鲁棒。
• 实操步骤（以DESeq2为例）：
  1. 用featureCounts（参数：-t exon -g gene_id -Q 30 --primary）生成raw count矩阵，确保只计数比对质量高（MAPQ≥30）、主比对（--primary）、外显子区域的reads；
  2. 导入R：dds <- DESeqDataSetFromMatrix(countData = counts_matrix, colData = sample_info, design = ~ condition)；
  3. 运行dds <- DESeq(dds)，DESeq2自动完成：a) 基于几何均值的size factor计算；b) 负二项模型拟合；c) Wald检验或LRT检验。
• 常见错误：把FPKM/TPM矩阵强行塞进DESeq2。我见过最离谱的案例：有人用TPM矩阵运行DESeq2，得到的log2FoldChange范围从-15到+20，而真实qPCR验证的最大变化只有±4倍。原因？TPM破坏了方差-均值关系，导致统计检验完全失效。

实操心得：featureCounts的-p（paired-end）和-B（require 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 '$3=="transcript" {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. 提取TPM（salmon输出的quant.sf文件第一列是Name，第四列是TPM） cut -f1,4 quant.sf > sample1.tpm

• FPKM/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直接告诉你“这个基因占转录组的百分之几”，数值越大，生物学意义越明确。例如TPM>100通常认为是高表达，TPM<1可能是低丰度或技术噪音。而raw_count受基因长度支配太大——一个100kb的胶原蛋白基因，raw_count=500可能只是基础表达；一个1kb的激酶基因，raw_count=500就是极高水平。
• 避坑指南：不要用TPM做“绝对定量”。TPM不是molecules/cell，它没有绝对物理单位。曾有个学生用TPM值去推算蛋白拷贝数，结果误差达3个数量级。TPM只适合相对比较（基因A vs 基因B，样本X vs 样本Y），不适合绝对丰度解读。

3.4 场景四：特殊实验类型——单细胞RNA-seq（scRNA-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供DEG，Salmon生成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 2：Hisat2比对（关键参数解析）

# 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

关键参数--dta（downstream transcript assembly）：告诉Hisat2保留所有比对信息，包括多比对位点，这对后续StringTie组装新转录本至关重要。不加此参数，StringTie会报错“no alignments found”。

Step 3：featureCounts生成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 'NR>1 {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 4：Salmon准确定量（转录本水平，生成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 # 提取TPM（quant.sf文件第四列） cut -f1,4 sample_salmon_quant/quant.sf > sample.tpm

实测对比：对同一组数据，featureCounts的基因计数与Salmon的TPM总和相关性达0.92，但Salmon在低丰度转录本（TPM<1）的检测灵敏度高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（用prcomp，scale. = 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 问题1：TPM矩阵中大量基因TPM=0，但raw_count显示有reads

现象：Salmon输出的quant.sf中，很多基因TPM=0，但featureCounts显示其raw_count>10。

根本原因：Salmon是转录本水平定量，TPM=0意味着没有足够证据支持该基因的任何转录本被表达。而featureCounts是基因水平计数，只要reads比对到该基因的任意外显子，就算入count。两者颗粒度不同。

排查步骤：

1. 用grep "ENSMUSG00000029197" Mus_musculus.GRCm39.104.gtf查看该基因的所有转录本ID；
2. 在Salmon的quant.sf中搜索这些转录本ID，看是否有非零TPM；
3. 若所有转录本TPM均为0，说明Salmon认为该基因无表达；若部分转录本TPM>0，则featureCounts的count可能来自未注释的转录本或比对错误。

解决方案：对于关注特定基因的研究，建议用featureCounts + StringTie联合流程：先用StringTie组装新转录本，再用featureCounts基于新GTF计数，最后用Ballgown做差异转录本分析。

5.2 问题2：FPKM/TPM值异常巨大（>10⁵）或为负数

现象：计算出的FPKM值高达500000，或TPM出现-0.0001。

根本原因：基因长度L_i为0或负数。常见于GTF文件中某些伪基因、lncRNA的exon坐标错误（start>end），或featureCounts的-g参数指定错误导致长度计算失败。

快速定位：

# 检查GTF中是否有start>end的行 awk '$4>$5' Mus_musculus.GRCm39.104.gtf | head -5 # 检查featureCounts输出的count文件，看是否有基因名为空或异常 head -10 sample.counts | cut -f1

修复方案：用gffread过滤GTF：

gffread -E Mus_musculus.GRCm39.104.gtf -o cleaned.gtf # -E 参数自动移除start>end的错误行

5.3 问题3：DESeq2运行报错“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) # 但此方法误差极大，仅用于快速预览，正式分析必须重跑featureCounts

5.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。