第4章 读取RNA-seq计数数据
开始做批次校正之前，得先把数据读进R。这一步听起来简单，但我在实际项目中见过太多人因为读取方式不对，后面分析全跑偏。你想想，一个基因表达矩阵如果基因ID对不上、列名错位、或者混进了非编码RNA，那后续所有结果都是垃圾。
我个人习惯先把原始计数数据整理成标准的矩阵或DGEList对象。edgeR的核心输入就是一个计数矩阵加上样本分组信息。今天我们就说说怎么把这堆txt、csv、或者featureCounts的输出，变成edgeR能认的格式。
4.1 常见计数数据格式
先搞清你手里的数据长什么样。下面三种最常见：
| 格式 | 来源 | 典型特征 |
|---|---|---|
| 基因×样本矩阵 | 自己整理 / 公共数据库 | 行是基因ID，列是样本名，数值是raw count |
| featureCounts输出 | Subread / Rsubread | 多列文件，包含基因ID、Chr、Start、End以及每个样本一列计数 |
| HTSeq-count输出 | HTSeq / STAR --quantMode | 每个样本一个文件，两列：基因ID和计数 |
大部分时候你会拿到featureCounts的表。嗯，这里要注意：表里可能有很多列，比如基因长度、外显子长度等，计数列只占其中一部分。读的时候不要一股脑全读进来，只留你需要的那几列。
4.2 用readDGE直接加载多个样本文件
如果你每个样本一个文件（比如HTSeq-count的输出），edgeR提供了readDGE函数。它可以把一堆文件合并成一个DGEList对象，非常方便。
我曾经接手一个项目，30个样本，每个样本一个单独的文件名还是乱的。我写了个循环把文件名整理好，一行readDGE搞定。避坑指南：文件里的基因ID列名必须统一，比如有的叫gene_id，有的叫GeneID……那就得先统一。
# 假设你的文件列表放在 files 向量里，每个文件两列：gene_id, count
library(edgeR)
files <- list.files(pattern = "*.counts$")
dge <- readDGE(files, columns = c(1,2))
# columns 指定基因ID列和计数列的位置（1-based）
默认readDGE会把文件名作为样本名。你可以用samples参数重新指定分组信息。
文件名最好包含分组和批次信息。例如
WT_batch1_sample1.counts，这样后面解析起来方便。我习惯用正则从文件名里提取分组和批次，省去手动输入。
4.3 手动构建DGEList对象
很多时候你拿到的是一个完整的计数矩阵（行基因，列样本）。那就不需要用readDGE了，直接自己建DGEList。
# 读取矩阵（假设是tab分隔）
raw_counts <- read.delim("count_matrix.txt", row.names = 1)
# 检查一下前几行
head(raw_counts)
# 构建DGEList
group <- factor(c("control","control","treatment","treatment"))
batch <- factor(c("batch1","batch2","batch1","batch2"))
dge <- DGEList(counts = raw_counts, group = group, genes = rownames(raw_counts))
# 把批次信息加到样本数据里
dge$samples$batch <- batch
为什么要手动加batch？因为edgeR的DGEList本身没有批次槽，我们后续建模时需要把它作为协变量放进设计矩阵。这个我们在第7章详细讲，但你先记住：批次信息必须和样本一一对应。
计数矩阵里不能有缺失值（NA），也不能有负数。如果你做的是双端测序，计数必须是整数。我曾经见过一个数据，里面混了FPKM……那根本不是raw count，千万别用。edgeR只接受整数raw count。
4.4 数据初步检查
数据读进来后，别急着往下跑。先做个快速检查：
- 总文库大小：用
colSums(dge$counts)看每个样本的总reads数，差距过大（比如10倍）说明可能有批次效应或样品问题。 - 基因检测数：每个样本里至少有一个count的基因数。如果某个样本基因数明显偏少，可能是建库失败。
- 基因ID去重：有些数据库会混入重复的基因名（比如ENSEMBL ID后面带了版本号）。我习惯用
rownames(dge$counts)检查重复值。
为什么要做这一步？因为一旦有错误数据，后面找差异基因时假阳性会飙升。我自己踩过坑：一个样本的计数单位弄错了（乘了100倍），结果差异表达列表里全是那个样本独有高表达的基因……浪费了一周排查。
下面这张图概括了从原始文件到DGEList的完整流程：
4.5 常见问题与解决
- 基因ID重复：用
duplicated()找出重复行，然后决定是取平均还是删掉。我个人倾向于保留第一个，但如果重复来自同一基因的不同转录本，建议合并计数。 - 计数列名包含特殊字符：R会自动把空格变成点，但最好在读取时用
check.names = FALSE，然后手动清理。 - 内存爆炸：RNA-seq数据动辄几万基因，几百样本时矩阵很大。用
data.table::fread()替代read.delim速度更快，内存更少。
读取数据不是简单的“读进来就行”。你需要确保原始计数是纯整数、基因ID唯一、并且能把分组和批次信息正确附加到DGEList中。这一步错了后面全白费。
好了，现在你手里已经有了一个干净的DGEList。下一件事就是检查批次效应是否存在。不过那是第5章的内容。先确认你的数据读对了吗？用dim(dge)看看维数，用dge$samples看看样本标签，确保一切就绪再往前推进。