📚 课程目录第2章安装与数据准备

2、安装与数据准备：导入count矩阵、colData，构建DESeqDataSet

好了，到了真正动手的环节。

DESeq2 的威力再大，也得先把数据喂进去才行。这一步说实话出错的概率不小——我自己刚入门时就因为 colData 的行名没对齐，排查了一下午。现在回想起来，就是少看了一眼文档。

这一章咱们把整个数据准备流程走一遍。从安装包开始，到构建出 DESeqDataSet 对象。搞明白每一步到底在做什么，后面调参数才心里有底。

2.1 安装 DESeq2 包

DESeq2 属于 Bioconductor 生态，不能用普通的 install.packages()。你得用 BiocManager 来装。

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("DESeq2")

注意，BiocManager 会自动帮你处理依赖关系。如果你的网络环境不太好，记得换个国内镜像。我个人习惯在 install.packages 时加个 repos 参数指向清华源，省心很多。

小提示： 装完后用 library(DESeq2) 测试一下，如果不报错，就说明安装成功。

2.2 准备 count 矩阵

DESeq2 要求的输入是原始整数计数，也就是每个基因在每个样本里测到的 reads 数。千万别用 RPKM / FPKM / TPM，那些已经做过长度和文库大小归一化了，DESeq2 内部自己会做标准化，你传进去反而坏事。

我曾在项目中接过一个合作方给的 FPKM 矩阵，完全没法跑。后来让他们重新导出了 featureCounts 的结果，才搞定。

count 矩阵的格式要求很简单：

行是基因（或转录本），列是样本。
每个值必须是正整数（≥0）。
列名必须和后面 colData 的行名一一对应。

一个典型的 count 矩阵长这样（假设是 CSV 文件读取后的样子）：

# 模拟的 count 矩阵（前 3 行，前 4 列）
         sample1 sample2 sample3 sample4
GeneA        152      89     230       45
GeneB          0      12       8        3
GeneC        678     432     901      210

实际项目中，count 矩阵往往来自 HTSeq-count、featureCounts 或 Salmon（用 tximport 转换）。我习惯把结果存成 counts.csv，然后用 read.csv() 读进来，把第一列设为行名。

countData <- read.csv("counts.csv", row.names = 1)

注意： 确保 countData 是 matrix 类型，而不是 data.frame。有些函数会隐式转换，但显式转一下更安全：countData <- as.matrix(countData)。

2.3 构建 colData —— 实验设计信息

colData 是个 data.frame，描述每个样本的实验分组和批次等信息。行名必须和 count 矩阵的列名完全一致（顺序可以不同，但名字必须一一对上）。

举个例子，如果我们要比较处理组和对照组：

# colData 示例
         condition  batch
sample1   control   batch1
sample2   control   batch1
sample3   treated   batch2
sample4   treated   batch2

创建一个 colData 的常规做法：

coldata <- data.frame(
    condition = factor(c("control", "control", "treated", "treated"),
                       levels = c("control", "treated")),
    batch = factor(c("batch1", "batch1", "batch2", "batch2")),
    row.names = c("sample1", "sample2", "sample3", "sample4")
)

这里我把 condition 转成了因子，并且明确设定了水平顺序。为什么要手动设？因为 DESeq2 默认会把第一个水平当作对照组（reference level）。不设的话，R 会按字母序排，万一 "treated" 排在 "control" 前面，结果就反了。我吃过这个亏，从此再也不敢偷懒。

还有一个容易忽略的点：colData 里除了分组信息，还可以加 batch、lane 等协变量。如果你的实验有明确的批次效应，一定要加进 colData，DESeq2 会在模型里自动校正。

2.4 构建 DESeqDataSet 对象

所有准备工作就绪后，用 DESeqDataSetFromMatrix() 函数整合起来。这一步相当于是把你的数据 + 实验设计 + 公式“打包”成一个专门的对象。

dds <- DESeqDataSetFromMatrix(
    countData = countData,
    colData = coldata,
    design = ~ batch + condition
)

design 公式指定了你要考虑哪些变量。格式和 R 的线性模型一样。如果你想控制批次，就加上 batch；如果只需要看分组差异，直接写 ~ condition 就行。

构建完成后，建议先检查一下对象的基本信息：

dds
# 会输出类似：
# class: DESeqDataSet 
# dim: 20000 4 
# metadata(1): version
# assays(1): counts
# rownames(20000): GeneA GeneB ...
# rowData names(0):
# colnames(4): sample1 sample2 sample3 sample4
# colData names(2): condition batch

看一眼assay和colData的内容是否和你预期一致。如果样本名没对上，会报错；如果没问题，就可以进入下一步分析了。

  关键检查点： 运行 rownames(colData(dds)) 和 colnames(counts(dds))，两者应该完全一致（顺序不重要，名字集合相同就行）。

构建好 dds 后，还有个实用操作：过滤低表达基因。我一般会先把 counts 矩阵里行总和小于 10 的基因去掉，能显著加快后续计算并减少多重检验校正的负担。

keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

过滤后 dds 的行数会减少很多，不过别担心，那些极低表达的基因本来就不可能检测到差异。

下面这张流程图展示了从原始数据到 DESeqDataSet 的完整路径：

看到这里，你已经把 DESeq2 的数据准备阶段摸透了。别小看这一步——我见过不少新手拿着漂亮的 count 矩阵，却因为 colData 写错一行，后面跑出来的结果完全没法解释。花点时间检查 all(rownames(colData(dds)) %in% colnames(counts(dds)))，值不值得？绝对值得。