📚 课程目录第2章数据准备与读取

数据准备与读取

做差异表达分析，第一步就是把数据喂给 edgeR。这一步看似简单，但我见过不少人在数据读入上翻车。说白了，数据准备没做好，后面分析全是白搭。这章我就带你走一遍数据准备与读取的完整流程，把坑都填上。

核心要点：edgeR 需要计数数据（整数）和分组信息。千万别给它标准化后的数据，也别给 RPKM/FPKM。我当年第一次跑 edgeR 就栽在这上面，教训深刻。

表达矩阵与分组信息

edgeR 接受的数据，本质上是一个基因×样本的计数矩阵。每一行是一个基因，每一列是一个样本。矩阵里的值必须是整数，代表比对到该基因的 read 数。

分组信息就是一个因子，告诉 edgeR 每个样本属于哪个处理组（比如对照组 vs 处理组）。这两个东西缺一不可。

我们来看一个典型的数据结构。假设你有6个样本，3个对照（Control_1, Control_2, Control_3），3个处理（Treatment_1, Treatment_2, Treatment_3）。表达矩阵大概长这样：

           Control_1 Control_2 Control_3 Treatment_1 Treatment_2 Treatment_3
GeneA          1520      1489      1632         1789        1812        1756
GeneB           234       210       198          567         589         543
GeneC            45        52        38          123         110         119
...

分组信息用向量表示：

group <- factor(c("Control","Control","Control",
                  "Treatment","Treatment","Treatment"))

嗯，这里要注意：因子的顺序决定了对照和处理组的基准。我习惯把对照组放在前面，这样默认比较就是处理 vs 对照。

常见数据格式

实际上我们拿到的数据很少是完美的矩阵。常见的来源有这些：

数据来源	典型格式	读入方式
HTSeq / featureCounts 输出	每个样本一个文件，两列：基因ID、count	readDGE() 合并
表达矩阵文本	CSV/TSV，基因行为名，样本为列	read.table() → DGEList()
DGEList 对象	R 对象，包含 counts、samples 等槽	直接加载 .rda 或 .rds
SummarizedExperiment	Bioconductor 标准容器	as.DGEList() 转换

你想想看，如果数据分散在几十个单独的文件里，手动合并多麻烦？edgeR 专门提供了 readDGE() 函数，一步到位。

使用 edgeR 读取数据

我们从最常见的场景开始。假设你已经用 HTSeq 跑完了定量，每个样本得到一个 _counts.txt 文件。文件名存在一个向量里：

files <- c("Control_1_counts.txt", "Control_2_counts.txt",
           "Control_3_counts.txt", "Treatment_1_counts.txt",
           "Treatment_2_counts.txt", "Treatment_3_counts.txt")
labels <- c("Control_1","Control_2","Control_3",
            "Treatment_1","Treatment_2","Treatment_3")

然后一行代码就能搞定：

library(edgeR)
x <- readDGE(files, labels = labels)
x$samples$group <- factor(c("Control","Control","Control",
                            "Treatment","Treatment","Treatment"))

我的小技巧：文件名字最好包含分组信息，比如用 list.files(pattern="*_counts.txt") 自动捕获。我习惯在文件名里加前缀（C_ 和 T_），这样分组一眼就能看出来。

如果数据已经合并成一个表达矩阵，比如叫 counts，那就直接用 DGEList()：

y <- DGEList(counts = counts, group = group)

注意 counts 必须是整数矩阵，不能有缺失值。我曾经遇到过一个数据集，里面居然有负数和 NA，读进来直接报错。这时候用 apply(counts, 2, as.integer) 强制转换，但更稳妥的是先排查数据来源。

⚠️ 常见坑：千万不要在 counts 矩阵中包含基因长度、GC含量等特征列。edgeR 只需要计数。如果你误把基因名当成一列，记得用 row.names() 设置行名后删掉那一列。

读入之后，DGEList 对象会自动生成一些统计量：

y$samples
   group lib.size norm.factors
Control_1 Control   1520000            1
Control_2 Control   1485000            1
Control_3 Control   1632000            1
Treatment_1 Treatment 1798000            1
Treatment_2 Treatment 1812000            1
Treatment_3 Treatment 1756000            1

lib.size 是每个样本的文库大小（总计数），norm.factors 初始为1，后面做 TMM 标准化时会更新。

为了让你对整体流程有个直观印象，我画了一张图：

两条路径最终都汇聚到 DGEList 对象。我个人更推荐文件列表的方式，因为能自动处理文件名和路径，不易出错。

检查读入的数据

数据读完了，别急着往下跑。先做几个快速检查：

维度：dim(y) 看看基因数和样本数对不对。
缺失值：any(is.na(y$counts)) 如果有 NA，必须处理。
非整数：all(y$counts == round(y$counts)) 确保没有小数。
分组平衡：table(y$samples$group) 检查各组样本数。

记住：如果 counts 矩阵里有小数，edgeR 会报错但不会告诉你具体原因。我遇到过一次，数据来自某个旧版定量工具，输出的是平均规范化后的计数。排查半天才发现问题——所以保险起见，每次读入后立刻检查数据类型。

好了，数据准备妥当，DGEList 对象已经就绪。下一阶段就是探索性分析和过滤低表达基因了——别着急，先把基础打牢，后面自然水到渠成。