📚 课程目录第3章数据输入格式（DESeq2）

3. 数据导入：读取计数矩阵和样本信息

3.1 数据从哪来？

做差异表达分析，第一步就是找对数据。说白了，你需要两样东西：计数矩阵和样本信息表。

计数矩阵长什么样？每一行是一个基因，每一列是一个样本，格子里的数字是比对到这个基因上的reads数。样本信息表呢？记录每个样本的分组、批次、测序lane这些实验条件。

我早期接过一个项目，对方甩给我一个Excel，基因名是中文，列名是日期……嗯，花了一整天清洗数据。所以我建议你从一开始就把格式固定下来。

数据类型	典型文件	行/列含义
计数矩阵	counts.csv, counts.txt	行=基因，列=样本
样本信息	colData.csv, metadata.txt	行=样本，列=分组/批次

3.2 用 read.table 之类的基础方法

最直接的做法：

counts <- read.table("counts.txt", header=TRUE, row.names=1)
colData <- read.table("metadata.txt", header=TRUE, row.names=1)

row.names=1 意思是用第一列作为行名。基因名通常放第一列。注意：CSV文件要用 read.csv，如果文件里有引号或逗号，记得加上 quote="\"。

警告： 我曾经在读取一个包含空格和特殊字符的基因名文件时，R 自动把空格变成了点。结果后续合并数据时处处报错。建议读完后马上检查 rownames(counts) 是否和你预期一致。

3.3 用 tximport 搞定定量结果

如果你用 Salmon、kallisto 或 RSEM 做的定量，那么你手里是转录本水平的估计，而不是基因水平。这时候别用 read.table 硬拼——用 tximport。

library(tximport)
files <- c("sample1/quant.sf", "sample2/quant.sf")
names(files) <- c("sample1", "sample2")
txi <- tximport(files, type="salmon", tx2gene=tx2gene, countsFromAbundance="no")

tx2gene 是一个从转录本ID到基因ID的映射表。为什么需要它？因为我们要把转录本水平的计数汇总到基因水平。我习惯在项目开始时用 ensembldb 或 biomaRt 造好这个映射，存成CSV，后面反复用。

小技巧： countsFromAbundance 参数有三个选项："no" 用原始估计计数，"scaledTPM" 或 "lengthScaledTPM" 做一些缩放。DESeq2 官方推荐用 "no"，因为它的内部自带校正。我试过 "scaledTPM"，结果和 "no" 的差异很小，但为了保险，还是用 "no" 吧。

3.4 样本信息表的艺术

你的 colData 必须要包含一个因子型的分组列。举个标准例子：

sample    condition    batch
sample1   treated      batch1
sample2   treated      batch1
sample3   control      batch2
sample4   control      batch2

确保行名和计数矩阵的列名完全一致。如果不一致，DESeq2会在 DESeqDataSetFromMatrix 时帮你对齐——但前提是两者都排序过，或者你用 match 手动调一下。我觉得“靠自己”更靠谱。

# 手动对齐样本顺序
counts <- counts[, rownames(colData)]

  核心原则： 计数矩阵的列名必须严格匹配 colData 的行名，而且顺序不重要，DESeq2 会按 colData 行名重新排列。

3.5 构建 DESeqDataSet 对象

这是数据导入的最后一步，也是最重要的桥梁。

library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData = colData,
                              design = ~ condition)

design 公式里的 ~ condition 表示你想比较不同condition组。如果有批次效应，要写成 ~ batch + condition。注意顺序——你关心的变量放最后。

同样，edgeR 也有对应的构建方式：

library(edgeR)
y <- DGEList(counts = counts, group = colData$condition)
y <- calcNormFactors(y)

两者的区别？DESeq2 把实验设计写在 design 里，edgeR 把分组信息写在 group 参数里。但核心思想一样：告诉软件“哪些样本是一组的”。

注意： edgeR 的 DGEList 不会自动检查样本名称匹配——你需要确保 counts 的列名和 group 的对应关系正确。我曾经把 sample1 和 sample2 的组别写反了，结果直接做差异分析，前十个基因全是假阳性。所以，读完后一定用 head(y$samples) 看一眼。

3.6 数据导入流程总览

3.7 常见坑与避坑指南

基因名重复：计数矩阵里不能有重复的行名。我之前在整合多个批次数据时，同一个基因被算了两次，导致 DESeq2 直接报错。用 duplicated(rownames(counts)) 检查。
非整数计数：DESeq2 只接受整数计数。如果 tximport 用了 countsFromAbundance="scaledTPM" 会得到小数，要小心。edgeR 可以接受小数，但经我测试最好还是用整数。
零计数过滤：很多人导入后直接跑差异分析，其实应该先过滤掉低表达基因。我会在构建对象前先过滤：counts <- counts[rowSums(counts >= 10) >= 3, ]，这样能减少多重检验的负担。
行名与列名自动转换：R 在读取 CSV 时可能把行名里的纯数字变成带 X 前缀的字符串。我习惯用 check.names=FALSE 避免。比如：

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

我的个人习惯： 所有原始数据文件都存放在 data/ 目录下，计数矩阵叫 counts_raw.csv，样本信息叫 colData.csv。写一个 R Markdown 脚本，从读取到构建对象一步完成。这样项目换人接手的时候，不用翻半天文件夹。

3.8 边缘情况：如何处理大矩阵？

假如你的基因数超过6万（比如鸡、猪的较新参考基因组），read.table 可能会很慢。我推荐用 data.table::fread 读，然后转成矩阵：

library(data.table)
counts <- as.matrix(fread("counts_big.csv"), rownames=1)

注意 fread 默认不会设置行名，你需要手动 rownames(counts) <- counts[,1] 再删掉第一列。或者用 fread(cmd="cut -f1 counts_big.csv") 先取行名。

另外，如果文件是 gzip 压缩的，fread 会自动处理，不用先解压。我遇到过磁盘空间不足，解压后文件太大，结果 fread 直接读压缩包完美解决。

好了，数据到手，接下来就是质量控制了。