📚 课程目录第4章 R语言绘图基础

4. 数据导入与探索：读取表达矩阵、样本元数据、基本统计量计算

数据导入是我每次做批次质控的第一步。说实话，这一步看着简单，但坑特别多。我刚入行时总以为read.table一把梭就够了，后来吃过不少苦头——文件编码不对、分隔符看走了眼、表头被当成了数据行……嗯，今天咱们把这些坑填上，再顺手把表达矩阵和元数据的基本统计量算清楚。

4.1 读取表达矩阵

表达矩阵通常是二维表，行是基因（或探针/蛋白），列是样本。常见的格式有CSV、TSV、TXT（制表符分隔）。我个人习惯用data.table::fread()，速度快，而且自动识别分隔符和表头，省心。

先看一个典型的TSV格式表达矩阵前几行长什么样：

GeneID  Sample1  Sample2  Sample3  Sample4
ENSG001  125.3    98.7     145.2    112.0
ENSG002  45.2     67.1     23.8     89.4
ENSG003  210.0    189.6    234.5    176.2
...       ...      ...      ...      ...

读取命令很简单：

library(data.table)
expr <- fread("expression_matrix.tsv")
dim(expr)
head(expr[, 1:6])   # 瞄一眼前几列

如果用的是CSV，记得指定sep = ","，但fread通常能自动嗅探出来。要是不放心，可以加上sep = "\t"显式声明。

小心行名被误读！ 很多表达矩阵第一列是基因名，但fread默认不设行名。你读进来后，第一列会是一个普通列。建议后续用expr <- as.data.frame(expr)然后手动设置rownames(expr) <- expr[[1]]；或者直接用read.table的row.names = 1参数。我两次因为忘了这步，后续画图时基因名全跑到了数据框里，排查了半天。

4.2 读取样本元数据

元数据（sample metadata）描述每个样本的批次、分组、处理条件等。它和表达矩阵是按列对应的——元数据的行是样本，表达矩阵的列也是样本，二者顺序必须一致，否则全乱套了。

我曾经接手一个项目，合作方发来的元数据和表达矩阵，样本排序完全相反，我用cbind合并后画出的PCA图，处理组和对照组完美分开——但其实是批次标签弄反了。所以，读取后第一件事就是检查样本名是否匹配。

元数据示例（CSV格式）：

SampleID,Batch,Treatment
Sample1,Batch1,Control
Sample2,Batch1,Treatment
Sample3,Batch2,Control
Sample4,Batch2,Treatment

读取：

meta <- fread("sample_metadata.csv")
meta <- as.data.frame(meta)
rownames(meta) <- meta$SampleID
meta$Batch <- factor(meta$Batch)
meta$Treatment <- factor(meta$Treatment)

这里我把Batch和Treatment转成了因子型——后续做分组统计或者绘图时，R会自动处理颜色和图例，省得我手动再转。

4.3 基本统计量——先摸清数据的脾气

拿到数据后别急着画图，先把几个关键数字撸出来。我一般会看这四样：

维度：多少基因 × 多少样本？
缺失情况：有没有NA、零值或者极低表达量？
表达量分布：中位数、四分位距、最大最小值。
样本间一致性：总表达量有没有某个样本特别低？

代码示例：

# 维度
cat("基因数:", nrow(expr), "样本数:", ncol(expr), "\n")

# 缺失值统计
na_count <- sum(is.na(expr))
cat("缺失值总数:", na_count, "\n")

# 每个样本的缺失值列数
na_per_sample <- colSums(is.na(expr))
cat("缺失样本列统计（前5）:\n")
print(head(sort(na_per_sample, decreasing = TRUE)))

# 表达量基本分布
library(matrixStats)
expr_mean <- rowMeans(expr, na.rm = TRUE)
expr_median <- rowMedians(as.matrix(expr), na.rm = TRUE)
quantile(expr_mean, probs = c(0, 0.25, 0.5, 0.75, 1))

至于样本间总表达量的比较，我通常做一个箱线图来快速筛查离群样本。不过在命令行阶段，可以先算每个样本的列和：

sample_total <- colSums(expr, na.rm = TRUE)
summary(sample_total)

一个小技巧：如果某个样本的总表达量比其它样本低一个数量级，别急着剔除。先检查原始数据——可能是这个样本的测序深度本身就低，也可能是批次效应。我遇到过类似情况，后来发现是文库制备时那个样本的扩增效率出了问题，属于技术批次因素，需要用后续的批次校正方法处理。

4.4 合并表达矩阵与元数据

做分析前，最好把两个数据结构的样本顺序对齐。我习惯用match()来保证列顺序与元数据行名一致：

# 确保表达矩阵的列名顺序与元数据行名相同
expr_ordered <- expr[, rownames(meta)]
# 检查是否全部匹配
all(colnames(expr_ordered) == rownames(meta))

如果返回FALSE，那就得找出不一致的样本，跟合作方确认。别偷懒跳过这步，我早期因跳过这步，后面画的热图标签全乱码了。

下面这张SVG图概括了本章的数据导入与探索流程，你跟着它做，基本不会漏掉关键节点。

4.5 实战：快速整理一个完整的探索报告

我把上面的步骤打包成一个函数，每次拿到新数据跑一遍就行。你会用到的代码：

explore_batch_data <- function(expr_path, meta_path) {
  # 读取
  expr <- fread(expr_path)
  meta <- fread(meta_path)
  # 转数据框并设行名
  expr <- as.data.frame(expr)
  rownames(expr) <- expr[[1]]
  expr[[1]] <- NULL
  meta <- as.data.frame(meta)
  rownames(meta) <- meta[[1]]
  meta[[1]] <- NULL
  # 对齐
  common <- intersect(colnames(expr), rownames(meta))
  expr <- expr[, common]
  meta <- meta[common, ]
  # 检查
  cat("表达矩阵:", nrow(expr), "x", ncol(expr), "\n")
  cat("元数据样本数:", nrow(meta), "\n")
  cat("缺失值总数:", sum(is.na(expr)), "\n")
  cat("表达量范围:", range(expr, na.rm = TRUE), "\n")
  cat("批次水平:", levels(factor(meta$Batch)), "\n")
  return(list(expr = expr, meta = meta))
}

把这个函数存放在你的R项目脚本库里，以后每次拿到新批次数据，一句res <- explore_batch_data("expr.tsv", "meta.csv")就搞定初步探索。

什么时候需要停下排查？ 如果发现以下问题，先别往下做可视化，回头检查原始数据：

缺失值占比超过5%（可以接受吗？要看检测平台和过滤策略）
某个样本的总表达量跟其他样本相差超过3倍标准差
批次变量里居然有未定义的标签（比如“Batch1”写成“ Batch1 ”带空格）
样本数对不上——表达矩阵有100列，元数据只有98行

别觉得这些琐碎。我记得有一次合作方给的元数据里，有一个样本ID末尾多了个回车符，R读进来后多了一个空行，导致match()返回NA。后来我用trimws()清洗列名才解决。所以读取之后，顺手做一下colnames(expr) <- trimws(colnames(expr))，能省去后续很多麻烦。

好了，数据导入和基本统计量这部分就到这里。你手头如果有自己的表达数据，建议现在就跑一遍上面的函数，看看输出的统计数字有没有异常。下一轮我们会把这些统计量画成图，直接判断批次效应的严重程度。