📚 课程目录第2章差异表达分析的基本流程

2 差异表达分析的基本流程

聊完工具选型的思路，咱们得先搞清楚差异表达分析到底在干嘛。说白了，就是从两组或多组样本里找出那些表达量有显著差异的基因。我刚开始接触这个流程的时候，总觉得“跑个软件就行”，后来踩了不少坑才明白，每个环节都可能让结果翻车。今天我把这套流程拆成几步讲清楚。

2.1 整体流程概览

差异表达分析不是一步到位的。它就像一条流水线，从原始数据进来，到最终拿到差异基因列表，中间得经过好几道工序。我习惯把这个流程画成下面这张图，你一看就明白。

好，我们一个一个环节拆开看。

2.2 数据预处理——别让垃圾数据毁了结果

原始FASTQ文件拿到手，第一件事不是直接跑定量，而是质量控制。我见过太多同学直接把低质量数据丢进去，结果差异基因列表里一堆假阳性。

我曾经吃过这个亏：有一批肿瘤样本，测序质量整体偏低，我贪省事没做trimming直接比对，最后差异分析跑出来一堆“显著差异”的基因，结果验证时发现全是测序错误的噪声。从此我再也不敢跳过质控。

质控常用的工具是FastQC和MultiQC。看完报告，重点关注：

每个碱基的测序质量（Q值）
接头序列污染
GC含量分布
重复序列比例

如果质量不好，就用Trimmomatic或cutadapt去接头、切低质量碱基。举个命令例子：

java -jar trimmomatic-0.39.jar PE \
  sample_R1.fastq sample_R2.fastq \
  sample_R1_clean.fastq sample_R1_unpaired.fastq \
  sample_R2_clean.fastq sample_R2_unpaired.fastq \
  ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \
  LEADING:20 TRAILING:20 SLIDINGWINDOW:4:20 MINLEN:36

质控完成后，下一步就是比对到参考基因组。常用的比对器有STAR、HISAT2、Salmon（比对模式是 quasi-mapping）。我个人习惯用STAR，速度够快，而且能直接输出每个基因的raw counts。

STAR --genomeDir /path/to/genome \
  --readFilesIn sample_R1_clean.fastq sample_R2_clean.fastq \
  --outSAMtype BAM SortedByCoordinate \
  --outWigType bedGraph \
  --outFileNamePrefix sample_

2.3 表达定量——从比对到计数矩阵

比对完得到BAM文件，需要用featureCounts或htseq-count来生成基因表达计数矩阵。featureCounts更快，而且是subread包的一部分，我一直用这个。

featureCounts -T 8 -p -t exon -g gene_id \
  -a annotation.gtf \
  -o counts.txt \
  sample1.bam sample2.bam sample3.bam

输出的counts.txt里每一列是一个样本，每一行是一个基因。如果想用转录本水平的定量，也可以用Salmon或Kallisto直接对FASTQ定量，省掉比对这一步——我个人的建议是：如果你只关心基因水平差异，Salmon足够；但如果想做可变剪接分析，最好还是STAR+featureCounts。

2.4 差异表达分析核心——DESeq2与edgeR

拿到计数矩阵之后，就到了最核心的一步：差异表达分析。目前最主流的工具是DESeq2（R包）和edgeR（R包）。它们都基于负二项分布模型，但细节上有点差别：

特性	DESeq2	edgeR
归一化方法	中位数比例	TMM
方差估计	经验贝叶斯收缩	经验贝叶斯+加权似然
支持小样本	强（尤其每群只有2-3个重复）	需要至少3个重复才稳定
输出结果	log2 fold change, pvalue, padj	logFC, PValue, FDR
典型代码复杂度	低（DESeqDataSetFromMatrix）	中（DGEList→calcNormFactors）

DESeq2的代码示例：

library(DESeq2)
countData <- read.table("counts_matrix.txt", header=TRUE, row.names=1)
colData <- read.table("sample_info.txt", header=TRUE, row.names=1)
dds <- DESeqDataSetFromMatrix(countData = countData,
                              colData = colData,
                              design = ~ condition)
dds <- DESeq(dds)
results <- results(dds, contrast=c("condition", "treat", "control"))
write.csv(as.data.frame(results), "DESeq2_results.csv")

edgeR的代码示例：

library(edgeR)
counts <- read.table("counts_matrix.txt", header=TRUE, row.names=1)
group <- c("control","control","treat","treat")
dge <- DGEList(counts=counts, group=group)
dge <- calcNormFactors(dge)
dge <- estimateDisp(dge)
fit <- glmQLFit(dge)
qlf <- glmQLFTest(fit, contrast=c(-1,1))
topTags(qlf, n=20)

一个小技巧：做差异分析前，先看看样本间聚类是否合理。我经常先用PCA或热图检查样本分组是否和预期一致。如果对照组和治疗组混在一起，那就算跑出差异基因也别太信。

2.5 多重检验校正——别被p值骗了

跑完差异性检验，你会得到每个基因的p值。但这里有个坑：人类基因组里大概有2万个表达基因，如果只拿p < 0.05筛选，纯随机也有1000个基因“显著”。所以必须做多重假设检验校正。

最常用的是Benjamini-Hochberg (BH) 方法，也就是控制假发现率（FDR）。DESeq2和edgeR默认都做这个校正，输出的padj或FDR就是校正后的p值。我一般用padj < 0.05且|log2FC| > 1作为筛选阈值。

我自己遇到过的一个典型情况：有一批数据，用原始p值筛选出800个差异基因，但校正后只剩下120个。老板一开始觉得“怎么这么少”，但验证实验表明120个里90%以上都是真的，而那些被砍掉的“假显著”基因里面有好多确实是噪声。

2.6 结果解读——最终筛选与可视化检查

得到差异基因列表后，我通常会做两件事：一是用散点图观察log2FC和p值分布，二是用火山图快速标记显著基因。火山图在后续章节会详细讲，这里先提一句，它横轴是log2 fold change，纵轴是-log10(padj)，非常直观。

除了火山图，我还喜欢用MA图（M=log2FC, A=average expression），用来检查差异分析是否受表达量水平影响。如果MA图上低表达基因的log2FC波动很大，说明你的归一化或滤波可能有问题。

最后，筛选差异基因时要注意：

不要只看padj，也要看log2FC的绝对值。有些基因padj很小但fold change只有1.1倍，生物学意义可能有限。
重复数量很重要。每组至少3个生物学重复，否则统计检验能力会很低。
记得做基因注释。得到基因ID后，用biomaRt或org.Hs.eg.db转成基因符号、名称。

好了，到此为止，差异表达分析的完整流程就串起来了。你想想看，从原始数据到差异基因列表，中间其实就这几步。但每一步的参数和细节都可能让结果天差地别。我建议你别着急，先把自己的数据跑一遍，每个环节都看看结果长什么样。