3. 构建 DGEList 对象
做差异表达分析，第一步就是把原始计数数据装进 edgeR 能认的格式。DGEList 就是这个格式的核心容器。说白了，它是一个打包好的数据对象，里面不光有表达量，还带着样本信息和基因注释。我刚开始用 edgeR 时，总觉得这步就是塞数据，后来发现其实不少坑藏在这里。
3.1 为什么需要 DGEList？
edgeR 内部的所有统计模型都建立在 DGEList 之上。它的设计思路很清晰——把数据、实验设计和基因属性绑在一起。你想想看，如果你每次都要把 count 矩阵、分组向量和基因长度分开传入，代码会变得多啰嗦？
- 统一管理：一个对象包含了 counts、samples、genes 三个核心组件。
- 自动对齐：当你想过滤低表达基因时，edgeR 会自动更新所有关联信息。
- 便于扩展：后续的 normalization、dispersion 估计都直接修改这个对象。
记得我去年帮一个同事调试 pipeline，他就是因为手动改了 count 矩阵但没同步更新样本分组，结果跑出来的差异基因全是错的。嗯，从那以后我坚持：只要数据有变化，就在 DGEList 里改，不要在外面另存一份。
3.2 DGEList 的核心结构
DGEList 本质上是一个包含多个数据框的列表。咱们拆开来看：
| 插槽（Slot） | 数据类型 | 作用 |
|---|---|---|
counts |
整数矩阵，行为基因，列为样本 | 原始 read/UMI 计数 |
samples |
数据框，行对应样本 | 样本分组、lib.size、norm.factors 等 |
genes |
数据框，行对应基因 | 基因 ID、符号、长度、GC 含量等 |
画个图看更直观：
这里有个关键点：counts 矩阵的行名必须与 genes 数据框的行名一致。我之前踩过坑，用 cbind 拼接时不小心把行名弄乱了，结果建出来的 DGEList 行数对不上，edgeR 直接报错。后来我学乖了，每次都用 identical(rownames(dge$counts), rownames(dge$genes)) 检查一遍。
3.3 从矩阵创建 DGEList
最常见的做法是从一个整数型的 count 矩阵开始。假设你已经有了一个矩阵 count_matrix：
library(edgeR)
# 最简单的情况 —— 只有 counts
dge <- DGEList(counts = count_matrix)
print(dge)
# 输出类似：
# An object of class "DGEList"
# $counts
# Sample1 Sample2 Sample3
# GeneA 10 15 12
# GeneB 0 3 1
# ...
# $samples
# group lib.size norm.factors
# Sample1 1 500000 1
# Sample2 1 620000 1
# Sample3 1 480000 1注意看 samples 数据框：默认 group 全是 1，lib.size 是自动算的列和，norm.factors 初始为 1。这时候你还没有指定分组，后续需要手动设置。
3.4 添加样本分组与基因注释
实际的实验设计往往有两组或更多组别。我一般推荐在创建 DGEList 时就传入分组信息，省得后面再改。做法很简单：
# 假设已有分组向量
group <- factor(c("Control", "Control", "Treatment", "Treatment"))
dge <- DGEList(
counts = count_matrix,
group = group
)
# 也可以创建后修改
dge$samples$group <- group基因注释可以用 genes 参数加入。比如你有一个 gene_annot 数据框：
# gene_annot 至少有 gene_id 和 symbol 两列
head(gene_annot)
# gene_id symbol length
# 1 ENSG001 TP53 2500
# 2 ENSG002 BRCA1 3200
dge <- DGEList(
counts = count_matrix,
group = group,
genes = gene_annot
)
# 验证一下
identical(rownames(dge$counts), rownames(dge$genes))
# [1] TRUEidentical 检查。另一个坑是基因符号中包含重复值，这会导致后续 filterByExpr 出错，最好提前用 make.unique 处理。
3.5 实战：构建一个完整的 DGEList
咱们用一个小型模拟数据走一遍完整的流程。这也是我平时测试新参数时的标准模板：
# 模拟 100 个基因，6 个样本
set.seed(42)
count_matrix <- matrix(
rnbinom(100*6, mu=100, size=5),
nrow = 100,
ncol = 6
)
rownames(count_matrix) <- paste0("Gene", 1:100)
colnames(count_matrix) <- paste0("Sample", 1:6)
# 分组：前3个对照，后3个处理
group <- factor(c(rep("Control", 3), rep("Treatment", 3)))
# 模拟基因注释
gene_annot <- data.frame(
Symbol = paste0("SYM", 1:100),
Length = round(runif(100, 1000, 5000)),
row.names = rownames(count_matrix)
)
# 构建 DGEList
dge <- DGEList(
counts = count_matrix,
group = group,
genes = gene_annot
)
# 查看结构
str(dge)
# List of 3
# $ counts : int [1:100, 1:6] ...
# $ samples:'data.frame': 6 obs. of 3 variables:
# ..$ group : Factor w/ 2 levels ...
# ..$ lib.size : num [1:6] ...
# ..$ norm.factors: num [1:6] 1 1 1 1 1 1
# $ genes :'data.frame': 100 obs. of 2 variables:
# ..$ Symbol: chr ...
# ..$ Length: num ...
# - attr(*, "class")= chr "DGEList"现在 dge 对象已经准备好，可以进入下一步的操作了。你可以用 dge$samples 看看 lib.size 是否合理：
dge$samples
# group lib.size norm.factors
# Sample1 Control 49528 1
# Sample2 Control 51203 1
# Sample3 Control 48765 1
# Sample4 Treatment 50891 1
# Sample5 Treatment 50134 1
# Sample6 Treatment 52010 1库大小都在 48k-52k 之间，分布均匀，没有明显异常。如果某个样本的 lib.size 明显偏小，比如只有别的样本的一半，那就要警惕了——可能是测序质量差或者计数矩阵有误。
3.6 常见陷阱与检查
- 计数矩阵非整数：RNA-seq 计数应该是整数（read count 或 UMI count）。如果你的矩阵里有小数，很可能是归一化后的值。edgeR 默认要求整数，但不会强制报错。我建议用
all(count_matrix == floor(count_matrix))先检查一下。 - 缺失值：edgeR 不支持 counts 里有
NA，必须提前处理。通常的做法是用0替代，或者直接删掉有 NA 的基因。 - 重复列名/行名：样本名和基因名必须唯一。重复名会导致 DGEList 内部索引混乱，edgeR 有时不报错，但结果会莫名其妙。用
anyDuplicated(rownames(count_matrix))检查。 - 分组顺序：默认情况下，
group因子的水平按字母或数字顺序排列。你需要注意哪个是对照组、哪个是处理组。我习惯显式指定参考组：group <- relevel(group, ref = "Control")。
检查清单：每次创建完 DGEList 后，运行下面三条命令，确认数据状态：
# 1. 检查 counts 有无 NA 或负数
summary(dge$counts)
# 2. 检查样本分组是否正确
table(dge$samples$group)
# 3. 检查基因注释行名是否匹配
head(rownames(dge$genes) == rownames(dge$counts))做到这步，你的 DGEList 就稳稳地建好了。接下来就可以用它做低表达基因过滤、计算 normalization factors，然后开始正式的差异分析。我个人觉得，把 DGEList 这块地基打扎实，后面 80% 的错误都能提前避免。你试试看？