2. 差异表达分析标准流程概述
开始讲差异表达分析之前，我想先聊聊整个流程。你可能在网上见过各种眼花缭乱的教程，说实话，核心骨架也就那么几块。我个人习惯把流程分成三大阶段：数据准备、定量比较、结果解读。每一步都有容易栽跟头的地方，我们挨个过一遍。
一句话总结：差异表达分析的目的是从两组或多组样品中找出表达水平有统计学显著差异的基因。标准流程就是一条从原始测序读到差异基因列表的流水线。
2.1 流程总览：从FASTQ到差异基因
下图是典型差异表达分析的标准流程。我标注了每个阶段的核心输入/输出，后面小节会逐一拆解。
2.2 每个阶段我踩过的坑
下面我结合自己跑过的项目，说说每个步骤最容易被忽略的细节。你回想一下，是不是经常拿到结果发现一堆假阳性？多半是某一步处理得不干净。
2.2.1 质控：别小看这一步
我见过有人跳过质控直接比对，结果比对率不到50%。说白了，低质量碱基和接头残留会拖垮后续所有步骤。我曾在一个RNA-seq项目中，样品测序质量分布尾部很差，用Trimmomatic切掉后reads损失了15%，但比对率从60%提升到92%。
我的习惯：每个样品单独跑FASTQC，把过表达序列、GC含量异常的样品单独标记出来。如果批次之间质控指标差异大，先排查实验原因，别急着往下走。
# 质控命令示例（使用Trimmomatic）
java -jar trimmomatic-0.39.jar PE \
sample_R1.fastq.gz sample_R2.fastq.gz \
sample_R1_trim.fastq.gz sample_R1_unpaired.fastq.gz \
sample_R2_trim.fastq.gz sample_R2_unpaired.fastq.gz \
ILLUMINACLIP:adapters.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:362.2.2 比对和定量：选对参考
基因表达定量现在主流有两种策略：比对后计数（如STAR+HTSeq）和轻量级比对（如Salmon/kallisto）。我个人的建议是：如果做差异表达，尽量用count-based方法（如STAR + featureCounts），因为大多数差异检测软件期望整数计数。曾经有位同学用Salmon的TPM直接跑DESeq2，结果发现DESeq2的DESeqDataSetFromMatrix要求输入原始计数，TPM是标准化后的值，这样会破坏统计模型。
| 定量工具 | 输入类型 | 适用场景 |
|---|---|---|
| STAR + HTSeq-count | FASTQ → BAM → count | 经典RNA-seq，标准差异分析 |
| Salmon / kallisto | FASTQ → transcript abundance (TPM) | 快速定量，需要后续换算为count |
| featureCounts | BAM → count | 替代HTSeq，运行更快 |
注意：不要直接把TPM或FPKM输入到DESeq2/edgeR中，它们的统计模型基于负二项分布，预期输入的是原始计数。如果用TPM，必须先用tximport或自己写的脚本转成近似count。
2.2.3 标准化：不是随便归一化就行
标准化这一步很多新手搞混。我们做差异分析用的是“缩放因子”标准化，目的是校正文库大小差异和组成效应。DESeq2默认的median-of-ratios方法比简单的CPM更稳健。我记得有一次项目里有两个样品的某几个高表达基因占了总count的40%，导致普通CPM把其他基因都压下去了。换成DESeq2的标准化后，差异基因列表完全变了。
# R语言示例：DESeq2标准化
library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = counts, colData = colData, design = ~ condition)
dds <- estimateSizeFactors(dds)
normalized_counts <- counts(dds, normalized=TRUE)2.2.4 差异检验：选择合适的方法
常用的工具有DESeq2、edgeR、limma（针对voom转化后的count）。如果你是小样本（n≤3），我建议用edgeR的exact test或DESeq2的Wald test，这两者都对小样本做了特殊处理。我曾经对比过DESeq2和limma在一个只有2个生物学重复的数据集上的表现：limma的voom对低表达基因的FDR控制偏松，而DESeq2相对保守一些。
2.3 避坑指南
做了几年差异分析，我总结了几条血泪教训：
- 批次效应不除，结果不可信。曾经有一次我合并了两批测序数据，直接用DESeq2跑差异，发现前50个差异基因全是批次标记基因。后来用
ComBat-seq校正后才得到真正的生物学差异。 - 不要只看p-value。我建议用log2 fold change（LFC）和adjusted p-value的双阈值，比如|LFC|>1，padj<0.05。纯靠p-value会筛出很多表达变化很微小但统计学显著的基因。
- 生物学重复至少3个。2个重复也能跑，但可靠性会大打折扣。我的学生曾经用2个重复跑出一堆差异基因，结果qPCR验证全部不显著。
- 注意低表达基因的过滤。通常在差异分析前应该过滤掉在所有样品中count<10的基因。不这么做会导致多重假设矫正的负担加重。
一句话：标准流程是死的，但数据质量是活的。每一步都做点检查，比如PCA看分组分离（下章会详细讲），这样跑出来的差异基因才经得起后续验证。
2.4 常用工具对比
下面这个表是我给团队培训时整理的，方便快速决策。
| 工具 | 统计模型 | 最适合场景 | 优点 | 缺点 |
|---|---|---|---|---|
| DESeq2 | 负二项，Wald/LRT | 多因素、时间序列 | 内置标准化、稳健 | 大样本时较慢 |
| edgeR | 负二项，exact test / GLM | 小样本、两组比较 | 速度快，灵活 | 需要手动处理设计矩阵 |
| limma + voom | 线性模型，加voom权重 | 微阵列、大样本RNA-seq | 处理复杂设计 | 对低计数敏感 |
| sleuth | 线性模型，Bootstrap不确定性 | 转录本级别，无比对定量 | 直接使用kallisto | 仅适用kallisto |
2.5 标准流程的GSEA延伸
跑完差异基因后，很多人直接跳到功能富集，但这里有个小陷阱：如果只用显著差异基因做GSEA（超几何检验），会丢失那些表达有趋势但未达阈值的基因。我建议保留全基因排序列表，用GSEA的preranked模式。这样信息量更大。
# 示例：用clusterProfiler做GSEA
library(clusterProfiler)
geneList <- sort(rank_vector, decreasing=TRUE)
gsea_result <- gseGO(geneList, OrgDb = org.Hs.eg.db, keyType = "ENTREZID", ont = "BP")