📚 课程目录第2章 DESeq2安装与数据准备（count矩阵、样本信息）

第二章：DESeq2安装与数据准备——别在起跑线上摔倒

安装DESeq2，说白了就一行代码的事。但我见过太多人在这个环节卡壳——不是因为不会装，而是装完之后数据格式不对，跑起来全是报错。

我个人习惯从一开始就把数据准备得整整齐齐，这样后面调参才会顺手。你想想看，如果连count矩阵都是歪的，再好的参数也救不了你对吧？

2.1 安装与加载

DESeq2托管在Bioconductor上，用官方的安装命令就好。注意，如果你是第一次装BiocManager，先把那个包装上。

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("DESeq2")

# 安装完成后加载
library(DESeq2)

一个小建议：装包时尽量用BiocManager::install()，它能自动处理依赖关系。别用install.packages("DESeq2")，那会从CRAN下旧版本，容易踩坑。

我在早期项目里就吃过这个亏——直接装了个CRAN上的0.x版，结果跟后来写的教程代码对不上，调试了整整一下午。后来学乖了，老老实实走Bioconductor渠道。

2.2 数据准备：count矩阵

DESeq2要求的count矩阵，核心就三点：

整数型——RNA-seq原始read counts，不是RPKM/FPKM，也不是log变换后的值。
行为基因/转录本，列为样本——这个顺序别搞反。
不能有缺失值——缺一个NA，DESeq2直接罢工。

举个实际例子，假设你有3个样本（WT1, WT2, KO1, KO2），count矩阵长这样：

# 一个示例count矩阵（前5行4列）
counts_matrix <- matrix(
  c(120, 98,  45,  66,
    230, 211, 312, 289,
    56,  43,  10,  8,
    889, 912, 754, 801,
    34,  27,  19,  23),
  nrow = 5,
  byrow = TRUE
)
rownames(counts_matrix) <- paste0("gene", 1:5)
colnames(counts_matrix) <- c("WT1", "WT2", "KO1", "KO2")
counts_matrix

输出类似这样：

      WT1 WT2 KO1 KO2
gene1 120 98   45  66
gene2 230 211 312 289
gene3 56  43   10   8
gene4 889 912 754 801
gene5 34  27   19  23

一个常见陷阱：如果你从featureCounts或HTSeq的输出读进来，往往是基因在列、样本在行。记得转置。我曾经在项目里忘了这一步，结果对每个样本跑出来的dispersion全是乱码，查了半天才发现是矩阵方向反了。

2.3 数据准备：样本信息表（colData）

样本信息表（也叫colData），行必须与count矩阵的列完全一致，列则记录实验设计的变量。比如处理条件、批次、配对关系等。

我们继续用上面4个样本的例子：

colData <- data.frame(
  row.names = c("WT1", "WT2", "KO1", "KO2"),
  condition = c("WT", "WT", "KO", "KO"),
  batch     = c("B1", "B2", "B1", "B2")
)
colData

输出：

     condition batch
WT1        WT    B1
WT2        WT    B2
KO1        KO    B1
KO2        KO    B2

这里有个细节：condition列最好是因子型，而且把对照组放在第一水平。当然DESeq2会自动处理，但手动设置因子顺序可以避免后续对比时拿错基线。

colData$condition <- factor(colData$condition, levels = c("WT", "KO"))

这样做的好处是，默认比较就是KO vs WT，符合大多数实验设计。

2.4 构建DESeqDataSet对象

数据和元数据准备好了，用DESeqDataSetFromMatrix()一把梭：

dds <- DESeqDataSetFromMatrix(
  countData = counts_matrix,
  colData   = colData,
  design    = ~ condition
)
dds

输出信息会显示class: DESeqDataSet，以及样本数、基因数、设计公式等。到这里，数据准备阶段就结束了，可以直接进入DESeq()标准流程。

核心流程示意：

2.5 数据准备的几个坑

空行或全零行：DESeq2会自动过滤掉所有样本里count总和小于1的基因，但你最好在dds <- dds[ rowSums(counts(dds)) > 1, ]之前自己检查一下。全零行留在里面虽然不影响结果，但会拖慢运行速度。
因子水平顺序：忘了设置因子顺序，DESeq2会按字母顺序把第一个作为基线。如果你的WT和KO刚好字母顺序相反，那结果对比方向就反了。我有一回给同事检查代码，发现他的log2FoldChange全部正负号反了，就是因为这个。
设计公式：最简单的~ condition只适合单因素实验。如果有批次效应，最好写成~ batch + condition。DESeq2会用线性模型把批次效应剥离掉，再找差异基因。
稀疏矩阵：DESeq2不支持稀疏矩阵格式（比如dgCMatrix）。如果你用Seurat或SingleCellExperiment导出数据，记得转成普通矩阵：as.matrix()。

一个小技巧：在正式跑DESeq()之前，先用assay(dds)看看count矩阵是否完整，用colData(dds)看看样本信息是否对齐。这样能省下后面debug的半小时。

2.6 本章小结

数据准备虽然看着简单，但往往是最容易出幺蛾子的环节。count矩阵和colData的行名列名一致性、因子水平顺序、全零行过滤——这三个点我建议你在每次跑分析前都过一遍。我自己写了个快速检查脚本，每次跑新数据前跑一遍，心里才踏实。