📚 课程目录第3章表达数据基础

表达数据基础：搞懂你的数据长什么样

先别急着跑模型，我们得先聊聊数据本身。说白了，limma 是个线性模型，它不关心你的样本是来自小鼠还是病人，它只认矩阵和注释。我见过太多人上来就做差异分析，结果一查表达矩阵，行名是乱的，列名是乱的，样本注释也对不上——那跑出来的结果你敢信吗？

这一节咱们就把数据这块地基打牢。你想想看，一个完整的表达数据对象，其实就三样东西：表达矩阵、样本注释、基因注释。搞清楚了这三样，后面的批次校正和差异分析才会顺风顺水。

表达矩阵：最核心的数据结构

表达矩阵说白了就是一个二维表格，行是基因/探针，列是样本。每个格子里的数值，可能是芯片信号强度，也可能是RNA-seq的count值或FPKM。我个人习惯把它想象成一个"基因在样本中的活性快照"。

先看一个典型的表达矩阵长什么样：

       样本1   样本2   样本3   样本4
基因A   12.5    10.2    15.8    9.1
基因B   5.3     6.7     4.2     7.9
基因C   20.1    18.4    22.6    19.8
基因D   8.9     11.3    7.6     10.5

注意几个硬性要求：

行名必须是基因ID或探针ID，不能有重复
列名必须是样本ID，要和样本注释文件里的ID一一对应
所有值必须是数值，不能有字符串或缺失

⚠️ 我在项目中遇到的坑： 有次合作方给的表达矩阵，样本名前面多了个看不见的空格，结果R读进来自动加了个"X"，跟样本注释文件对不上，白白排查了半天。所以读完数据后，第一件事就是用 colnames() 和 rownames() 检查一下。

样本注释文件：给样本贴上标签

表达矩阵只管数值，样本注释文件则告诉我们每个样本的"身份信息"——是处理组还是对照组？是批次1还是批次2？这玩意儿在limma里叫 phenoData，其实就是个数据框。

一个标准的样本注释文件大概长这样：

sampleID    group    batch    sex
样本1      control  1        M
样本2      control  1        F
样本3      treat    2        M
样本4      treat    2        F

这里有几个关键点：

第一列必须是样本ID，而且要和表达矩阵的列名完全一致
分组变量要声明成因子，R才不会自作主张按字母顺序排序
批次变量也建议转成因子，虽然数值也能用，但因子更安全

读取数据：read.table 的实战技巧

真实场景下，表达矩阵和样本注释通常分开存储在两个文件里。读取它们的方式很简单，但有讲究。我个人偏好的做法是：

# 读取表达矩阵
expr_matrix <- read.table("表达矩阵.txt", 
                        header = TRUE, 
                        row.names = 1,
                        sep = "\t",
                        check.names = FALSE)

# 读取样本注释
pheno_data <- read.table("样本注释.txt", 
                       header = TRUE,
                       row.names = 1,
                       sep = "\t",
                       stringsAsFactors = FALSE)

💡 三个关键参数：

header = TRUE：把第一行当列名
row.names = 1：把第一列当行名
check.names = FALSE：防止R把样本名中的横线或空格改成点号

读完之后一定要做的事情——验证数据一致性：

# 检查样本ID是否一致
all(colnames(expr_matrix) == rownames(pheno_data))

# 如果不一致，找出差异
setdiff(colnames(expr_matrix), rownames(pheno_data))

构建 ExpressionSet：把数据整理成标准格式

在Bioconductor的世界里，ExpressionSet 是表达数据的标准容器。它把表达矩阵、样本注释、基因注释打包成一个对象，方便后续操作。我刚开始用Bioconductor时觉得这玩意儿多此一举，后来发现它的好处——你不需要记住三个数据框的名字，一个对象全搞定。

构建 ExpressionSet 的完整流程：

library(Biobase)

# 1. 创建样本注释对象 - AnnotatedDataFrame
pheno_df <- new("AnnotatedDataFrame", data = pheno_data)

# 2. 创建基因注释 - 可选，但推荐
# 假设我们有一个基因注释数据框 gene_anno
# 行名是基因ID，列是基因符号、染色体位置等
feature_df <- new("AnnotatedDataFrame", data = gene_anno)

# 3. 创建 ExpressionSet
eset <- ExpressionSet(
    assayData = as.matrix(expr_matrix),
    phenoData = pheno_df,
    featureData = feature_df
)

# 查看对象
eset

构建完成后，你会发现所有信息都在一个对象里：

exprs(eset) —— 取表达矩阵
pData(eset) —— 取样本注释
fData(eset) —— 取基因注释

组件	访问函数	说明
表达矩阵	`exprs(eset)`	行=基因，列=样本
样本注释	`pData(eset)`	行=样本ID，列=分组/批次等变量
基因注释	`fData(eset)`	行=基因ID，列=基因符号/位置等
实验描述	`experimentData(eset)`	实验设计、引用等信息

在构建过程中还遇到过哪些坑

我调一个真实案例——有次从GEO下载表达数据，样本注释文件的样本名排序跟表达矩阵不一样。你猜怎么着？直接构建ExpressionSet也没报错，但后续分析就全乱了。为什么？因为ExpressionSet在构建时会自动按样本名排序，但如果文件名里有特殊字符，排序结果可能和你想的不一样。

所以我在构建前必做的一步是：

# 确保两个文件的行名（样本名）排序一致
pheno_data <- pheno_data[colnames(expr_matrix), , drop = FALSE]

# 然后再构建
pheno_df <- new("AnnotatedDataFrame", data = pheno_data)
eset <- ExpressionSet(assayData = as.matrix(expr_matrix), 
                      phenoData = pheno_df)

✅ 安全构建三步走：

表达矩阵的列名 = 样本注释的行名
样本注释按表达矩阵列名排序
转成矩阵再构建（as.matrix）

实战经验小结

说了这么多，其实就三件事：

数据读进来要验证，别信文件里的数据一定对
样本ID是纽带，表达矩阵和注释文件靠它连起来
构建ExpressionSet是一劳永逸的事情，后面做分析直接调对象，不用再翻原始文件

我个人习惯于把读数据和构建这步写成函数，每次新项目直接调用。代码重复利用，也减少出错机会。嗯，准备好这些，我们才能放心地把数据交给limma去做批次校正。