1. RNA-seq 差异分析简介
1.1 什么是 RNA-seq 差异分析
RNA-seq 差异分析，说白了就是比较不同条件下基因表达到底变了多少。你拿到一批测序数据，想知道处理组和对照组之间哪些基因上调了、哪些下调了，这就是差异分析要回答的核心问题。
我 刚开始接触RNA-seq时，总觉得这步就是个统计检验，跑个软件就行。后来才发现，理解背后的逻辑比敲命令更重要。你想想看，一个基因在样本A里表达量是100，在样本B里是200，这算不算差异？如果A和B各只有两个重复，这个差异到底靠不靠谱？这些都需要差异分析来回答。
差异分析的应用场景其实很广：
- 疾病 vs 正常：找致病基因或标志物
- 给药 vs 对照：评估药物靶点和药效
- 不同发育阶段：研究调控机制
- 敲除 vs 野生型：验证基因功能
核心概念：差异表达分析（Differential Expression Analysis）是通过统计假设检验，识别在不同实验条件下表达水平有显著差异的基因。它的输出通常是一个基因列表，包含 fold change、p-value、adjusted p-value 等关键指标。
你可能会问，为什么不直接用倍数变化来筛选？嗯，这里有个坑。我 以前做过一个项目，用纯 fold change 筛基因，结果挑出来的全是低表达基因，噪声极大。后来才明白，必须结合统计检验，用 p-value 和 FDR 来控制假阳性。
1.2 差异分析的核心工作流
差异分析不是从原始数据直接跳到基因列表。它背后有一套完整的工作流。我 个人习惯 把流程分成六个阶段，每一步都影响最终结果的质量。
每个阶段都有自己的关键任务：
- 原始数据：FASTQ 格式，包含 reads 序列和质量分数
- 质控：检查数据质量，过滤低质量 reads，去除接头
- 比对：将 reads 比对到参考基因组或转录组
- 定量：统计每个基因的 read 计数，生成表达矩阵
- 差异分析：使用统计模型识别差异表达基因
- 可视化：用图形展示结果，辅助解读
我的经验：质控这一步千万别省。我 曾经遇到过 一批数据，差异分析结果很奇怪，折腾了两天才发现是接头污染导致的。后来我养成了习惯，每次拿到数据先跑 FastQC，再跑 MultiQC 汇总，一目了然。
1.3 常用差异分析工具
目前主流的差异分析工具就那么几个，但各有各的脾气。我先列个表，再一个个聊。
| 工具 | 方法原理 | 适用场景 | 优点 | 注意点 |
|---|---|---|---|---|
| DESeq2 | 负二项分布 + 收缩估计 | 小样本、多重复 | 稳健、假阳性控制好 | 对异常值敏感 |
| edgeR | 负二项分布 + 经验贝叶斯 | 中等样本量 | 速度快、灵活 | 需手动归一化 |
| limma | 线性模型 + 经验贝叶斯 | 大样本、微阵列出身 | 计算快、易扩展 | 对 count 数据需转换 |
| DEXSeq | 外显子级别模型 | 差异外显子使用 | 精度高 | 计算耗时长 |
说实话，我最常用的是 DESeq2。为什么？因为它对新手最友好。归一化、差异检验、多重检验校正一条龙，代码量少，出图也方便。我 用过的工具中，DESeq2 的文档也是最清楚的。
但 edgeR 也有它的优势。比如你要自定义对比矩阵、做复杂的实验设计，edgeR 的灵活性就体现出来了。我 记得有个项目 需要做 time-course 分析，edgeR 的 glm 框架用起来很顺手。
至于 limma，它本来是为芯片数据设计的，后来也支持 RNA-seq。它最大的优点是快，适合样本量大的场景。但你得先用 voom 把 count 数据转成 logCPM，这一步很多人会忽略。
来看一个 DESeq2 的基本用法：
library(DESeq2)
# 构建 DESeq2 对象
dds <- DESeqDataSetFromMatrix(
countData = count_matrix,
colData = col_data,
design = ~ condition
)
# 运行差异分析
dds <- DESeq(dds)
# 提取结果
res <- results(dds, contrast = c("condition", "treated", "control"))
res <- as.data.frame(res)
# 筛选差异基因
deg <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
选型建议：我个人习惯是，样本量少（< 8 个）用 DESeq2，样本量大（> 20 个）用 limma，需要灵活建模用 edgeR。当然，这不是绝对的，还得结合数据特点来定。
1.4 可视化——让差异"看得见"
差异分析跑完后，你会得到一个基因列表，里面有上千个基因的 p-value 和 fold change。但光看表格，你很难快速判断数据质量、整体分布和生物学意义。这时候就需要可视化登场了。
我 在项目中经常遇到 这样的情况：统计检验结果很显著，但一看可视化图，发现数据分布有问题或者存在批次效应。所以可视化不只是"画图交差"，它是分析质量的"照妖镜"。
最常见的几种差异可视化手段：
- MA 图：展示表达量（A）与倍数变化（M）的关系，用来判断数据整体分布
- 火山图（Volcano plot）：结合 fold change 和 p-value，直观展示差异基因
- 热图（Heatmap）：展示差异基因的表达模式，做样本聚类
- PCA 图：在差异分析前使用，评估样本分组和批次效应
举个简单的例子，画一张火山图：
library(ggplot2)
# 准备数据
res$gene <- rownames(res)
res$significant <- "No"
res$significant[res$padj < 0.05 & res$log2FoldChange > 1] <- "Up"
res$significant[res$padj < 0.05 & res$log2FoldChange < -1] <- "Down"
# 绘图
ggplot(res, aes(x = log2FoldChange, y = -log10(padj), color = significant)) +
geom_point(alpha = 0.6, size = 1.5) +
scale_color_manual(values = c("Up" = "#e74c3c", "Down" = "#3498db", "No" = "#bdc3c7")) +
theme_minimal() +
labs(title = "Volcano Plot", x = "log2(Fold Change)", y = "-log10(adjusted p-value)")
避坑提醒：我 曾经踩过的一个坑：用 p-value 而不是 adjusted p-value 来画火山图。结果出来一堆假阳性，差点把无关基因当成重要靶点。后来导师一眼就看出来了，说"你这个 p-value 没校正吧？" 场面一度很尴尬。所以，一定用 padj 或 FDR，别偷懒。
可视化的真正价值在于：它能把枯燥的数字变成直观的模式。好的可视化能让你一眼看出数据质量好不好、分组是否合理、差异基因是否集中在某些通路。这些信息，光看表格是看不出来的。
所以，差异分析和可视化是"一条龙"的关系。分析结果需要可视化来验证和解读，而可视化又能反过来指导分析策略的调整。两者缺一不可。