📚 课程目录第3章转录组学与基因表达图谱

3. 转录组学与基因表达图谱

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值分布直方图。如果分布不对劲，那就停下排查问题，别往下跑了。