第4章：R基础回顾：数据框、矩阵操作、读取单细胞数据

在正式使用 Harmony 校正批次效应之前，我建议大家先把手底下的 R 基本功夯实一下。

说实话，单细胞数据分析中绕不开的就是数据框和矩阵的来回倒腾。稀疏矩阵、表达矩阵、meta.data，这些结构你天天都要打交道。如果基础不熟，后面 Harmony 跑完了你都不知道结果是不是对的。

我最早做单细胞项目时，就因为搞混了data.frame和matrix的索引方式，结果抽了一个小时的细胞亚群，全抽错了……后来我就养成了一个习惯：每次拿到新数据，先检查数据类型和维度，再往下走。

今天咱们就快速过一遍数据框、矩阵的核心操作，再演示如何把 10X 的数据读进 R 里。文末还有一张数据处理流程图，建议你先看一遍再动手。

4.1 数据框操作——你每天都要用的那几个

数据框在单细胞分析里最常见的地方就是meta.data，也就是细胞的注释表。说实话，大部分操作就是增删改查。

4.1.1 创建与查看

# 创建一个简单的 meta.data 数据框
meta_df <- data.frame(
  cell_barcode = paste0("cell_", 1:10),
  sample_id    = rep(c("sample_A", "sample_B"), each = 5),
  nCount_RNA   = round(runif(10, 1000, 5000)),
  percent.mt   = round(runif(10, 2, 15), 1),
  stringsAsFactors = FALSE
)

# 看一眼
head(meta_df)
str(meta_df)

4.1.2 选取行/列

数据框取子集，我个人偏爱用$和[ , ]的组合，不过如果你追求速度，可以学学dplyr的filter和select。

# 选取列
meta_df$sample_id                      # 返回向量
meta_df[, "nCount_RNA"]                # 返回向量
meta_df[, c("cell_barcode", "sample_id")]  # 返回数据框

# 选取行（条件过滤）
high_umi <- meta_df[meta_df$nCount_RNA > 3000, ]
high_umi

4.1.3 添加新列

# 给数据框加一列，比如把样本名改为大写
meta_df$SAMPLE <- toupper(meta_df$sample_id)
head(meta_df)

4.1.4 合并数据框

项目里经常需要把多个样本的 meta.data 拼在一起。横向合并用merge，纵向合并用rbind。

# 纵向合并（要求列名完全一致）
meta_merged <- rbind(meta_df, meta_df)
dim(meta_merged)

# 横向合并（以 cell_barcode 为键）
extra_info <- data.frame(
  cell_barcode = meta_df$cell_barcode,
  cluster_id   = sample(1:5, 10, replace = TRUE)
)
meta_merged2 <- merge(meta_df, extra_info, by = "cell_barcode")
head(meta_merged2)

4.2 矩阵操作——单细胞表达数据的核心

单细胞表达矩阵通常是稀疏矩阵（dgCMatrix），但为了讲清楚基础，咱们先用常规矩阵演示。记住：行是基因，列是细胞，这是 Seurat 的惯例。

4.2.1 创建矩阵

# 模拟一个小表达矩阵：100 个基因 × 50 个细胞
set.seed(42)
expr_mat <- matrix(
  round(runif(100*50, 0, 10)),
  nrow = 100, ncol = 50,
  dimnames = list(
    paste0("Gene", 1:100),
    paste0("Cell", 1:50)
  )
)
expr_mat[1:5, 1:5]

4.2.2 维度与索引

dim(expr_mat)      # 100 50
nrow(expr_mat)     # 100
ncol(expr_mat)     # 50

# 取前10个基因，前5个细胞
small_mat <- expr_mat[1:10, 1:5]

# 取某个基因在所有细胞的表达
expr_mat["Gene42", ]

# 取某个细胞的所有基因表达
expr_mat[, "Cell13"]

4.2.3 矩阵转置

为什么要提转置？因为Harmony 输入的数据格式是 cells × PCs，而 Seurat 里的降维结果通常是 rows = cells，不需要转置；但如果你自己算 PCA，就要小心方向。

# 转置后变成：50 个细胞 × 100 个基因
t_expr <- t(expr_mat)
dim(t_expr)

4.2.4 矩阵 vs 数据框：对比清单

4.3 读取单细胞数据——坐上驾驶位

读取 10X 数据是每个单细胞分析师的必修课。这里我用Seurat::Read10X()演示，其实底层就是读三个文件：matrix.mtx、barcodes.tsv、features.tsv。它们组合起来就是一个稀疏矩阵。

# 假设你已经下载了 10X 样本数据到 ./data/sample1/
library(Seurat)

raw_counts <- Read10X(data.dir = "./data/sample1/")
class(raw_counts)   # dgCMatrix
dim(raw_counts)     # 基因数 × 细胞数

# 快速查看前几个基因在头几个细胞中的表达
raw_counts[1:5, 1:5]

4.3.1 从矩阵到 Seurat 对象

读进来的稀疏矩阵，直接喂给CreateSeuratObject即可。

seurat_obj <- CreateSeuratObject(
  counts = raw_counts,
  project = "sample1",
  min.cells = 3,     # 至少在3个细胞中表达的基因才保留
  min.features = 200 # 每个细胞至少检测到200个基因
)

seurat_obj
# 查看 meta.data
head(seurat_obj@meta.data)

4.3.2 手动构造 —— 为了理解数据结构

我建议新手至少手动造一遍 Seurat 对象，这样后面改数据不会心虚。

# 手动创建一个迷你 Seurat 对象
mini_counts <- raw_counts[1:100, 1:20]
mini_meta <- data.frame(
  row.names = colnames(mini_counts),
  sample_id = "sample1",
  batch     = "A"
)
mini_obj <- CreateSeuratObject(
  counts = mini_counts,
  meta.data = mini_meta
)
mini_obj

4.4 本章知识体系一览

4.5 实战避坑清单（我踩过的坑）

• 行/列方向搞反：Harmony 默认输入是 cells × PCs，而 Seurat 的降维矩阵是 cells × dimensions，方向天然对齐；但如果你自己用 prcomp() 算 PCA，记得转置。我曾经直接拿基因表达矩阵（基因×细胞）输进去，结果 Harmony 跑出负数……
• meta.data 中的批次列名不要带空格或特殊符号：Harmony 内部会按你指定的列名分组，如果列名里有反引号，容易报语法错误。我现在的习惯：一律用下划线，比如 batch_id。
• 稀疏矩阵转 data.frame 时内存爆炸：如果你不小心用 as.data.frame(expr_mat) 把整个表达矩阵转成了数据框，10 万细胞 × 2 万基因的数据会让你 R session 直接挂掉。用 as.matrix() 也会崩。正确的做法：按需取子集，或者用 tidyverse 的 as_tibble() 配合 rownames_to_column() 在稀疏矩阵上操作。
• 使用 merge 合并 meta.data 时行顺序被打乱：merge 默认会按 key 排序，导致合并后的数据框行顺序与原始 Seurat 对象细胞顺序不一致。所以我推荐用 dplyr::left_join 或者手动按行名匹配，总之要保证 row.names 对齐。