第3章 转录组学与基因表达图谱：RNA-seq技术、差异表达分析与可视化
3.1 RNA-seq技术：从样本到计数矩阵
说实话，RNA-seq现在几乎成了每个转录组项目的标配。它不像芯片那样需要预先设计探针，而是直接对整个转录组进行测序。说白了，就是把细胞里所有的RNA分子反转录成cDNA，然后上机测序，最后通过比对和定量，得到每个基因的表达量。
我个人习惯把这个流程拆成三步：
- 建库与测序：mRNA富集或rRNA去除，然后片段化、反转录、加接头。
- 比对与定量：用STAR或HISAT2将reads比对到参考基因组，再用featureCounts或HTSeq-count得到基因层面的counts矩阵。
- 质量控制：比对率、基因体覆盖均匀度、rRNA污染比例——这些指标一个都不能少。
避坑指南：我曾经接手过一个RNA-seq项目，测序深度高达100M reads，但比对率只有40%。后来发现是样本提取时混入了大量gDNA。所以拿到fastq文件后，第一件事是用FastQC检查GC含量和duplication rate，别急着往下跑。
常用的比对工具我就列一个简单对比：
| 工具 | 速度 | 内存需求 | 适合场景 |
|---|---|---|---|
| STAR | 极快 | ~30GB | 标准RNA-seq，高精度 |
| HISAT2 | 快 | ~5GB | 资源受限环境，大型项目 |
| Salmon | 极快 | ~8GB | 无参考基因组时的定量 |
注意：多组学项目中，转录组的数据量往往是最大的。一个常规的人类转录组项目，原始fastq文件动辄几十GB。建议提前规划好存储和计算资源，不然你想想看——跑到一半磁盘满了，那才真是欲哭无泪。
3.2 差异表达分析：找到真正的信号
拿到counts矩阵后，下一步就是找差异表达基因。这里有个问题——基因表达数据是离散的计数，服从负二项分布，所以不能用普通的t检验。
我在项目中用得最多的是DESeq2。它对小样本量更稳健，内置了离群值检测和多重检验校正。另一个常用工具是edgeR，它的速度更快，适合大样本量数据集。
差异表达分析的关键前提：
- 每个基因至少有10个以上的总counts。
- 至少3个生物学重复，否则不能估计离散度。
- 批次效应必须提前校正，否则结果根本不可信。
下面是DESeq2的标准用法：
library(DESeq2)
# 构建DESeqDataSet
dds <- DESeqDataSetFromMatrix(
countData = count_matrix,
colData = col_data,
design = ~ condition
)
# 运行差异表达分析
dds <- DESeq(dds)
# 提取结果
res <- results(dds, contrast = c("condition", "treatment", "control"))
res <- res[order(res$padj), ]
嗯，这里要注意——log2 fold change的收缩（shrinkage）很重要。我见过很多新手直接用原始log2FC排序，结果在低表达基因上选出大量假阳性。用lfcShrink函数处理后，结果会稳定很多。
3.3 热图与火山图：用眼睛找模式
差异表达分析结束后，通常会得到几百到几千个候选基因。这时候，可视化就派上用场了。我常用的两个图是热图和火山图。
热图用于展示样本间的表达模式。行是基因，列是样本，颜色深浅代表表达量高低。我一般会用pheatmap包来画。
# 提取显著差异基因的表达矩阵
sig_genes <- rownames(res[res$padj < 0.05, ])
mat <- assay(rld)[sig_genes, ] # rld是经过rlog变换后的数据
# 绘制热图
pheatmap(
mat,
scale = "row",
annotation_col = col_data[, "condition", drop = FALSE],
show_rownames = FALSE,
fontsize_row = 6
)
火山图则更直观地展示差异表达的整体情况。横轴是log2FC，纵轴是-log10(padj)。左上角是下调基因，右上角是上调基因。我用EnhancedVolcano包快速出图：
library(EnhancedVolcano)
EnhancedVolcano(
res,
lab = rownames(res),
x = "log2FoldChange",
y = "padj",
pCutoff = 0.05,
FCcutoff = 1,
colAlpha = 0.5
)
我的一点经验：火山图有个坑——当样本量小时，padj分布会不够均匀，导致图上一堆点挤在底部。我建议先做一个MA图看看，如果MA图的散点分布不对称，那可能是批次效应没处理好，不要急着下结论。
3.4 本章知识体系
我把转录组图谱分析的核心流程总结成一张图。你从左上角开始，顺着箭头走一遍，就知道每一步该做什么。
这张图里，从样本准备到生物学解释一共五步。每一步的输出都是下一步的输入。你想想看，如果中间任何一步出了问题，最终结果都可能跑偏。
3.5 常见问题与避坑
最后我总结几个实际项目中踩过的坑：
- 样本数量不足：别信什么2个重复就能做差异分析。我测试过，2 vs 2的p值分布是均匀的，基本没有统计学效力。至少3 vs 3起步。
- 忽略批次效应：如果你的样本是在不同时间做的文库，那批次效应几乎必然存在。用sva包的ComBat函数或者limma的removeBatchEffect。
- 只看p值不看效应量：一个基因p值很小但log2FC只有0.1，这在生物学上基本没意义。我个人习惯先设一个|log2FC| > 1的阈值。
- 热图颜色尺度不对：如果不做行缩放（scale="row"），高表达基因会把颜色拉偏，低表达基因的差异就看不出来了。
最后再提醒一次：转录组数据量大、分析步骤多，每一步的结果都可能影响下游。建议每做完一步就检查一下输出。我自己的习惯是——比对完立刻看比对率，定量完立刻看基因覆盖度，差异分析完立刻看p值分布直方图。如果分布不对劲，那就停下排查问题，别往下跑了。