📚 课程目录第2章数据准备

数据准备：读取10X Genomics数据、创建Seurat对象、数据格式检查

做单细胞分析的第一步，就是把原始数据请进R环境里。这一步要是没走稳，后面什么降维聚类都是白搭。我见过不少初学者一上来就直接跑NormalizeData，结果报错半天才发现是数据没读对。今天咱们就稳扎稳打，把读取10X数据、创建Seurat对象、检查格式这三件事一次说透。

2.1 读取10X Genomics数据

10X的数据输出通常是三个核心文件：matrix.mtx（稀疏矩阵）、barcodes.tsv（细胞条形码）、features.tsv（基因名）。它们都放在一个叫filtered_feature_bc_matrix（或raw_feature_bc_matrix）的文件夹里。我个人习惯把所有样本的原始数据按样本名/outs/filtered_feature_bc_matrix/的路径存放，这样后期批量处理时非常方便。

Seurat专门为10X数据写了个封装函数——Read10X。它自动识别三个文件并返回一个稀疏矩阵。来，看代码：

library(Seurat)

# 指定10X输出目录
data_dir <- "data/sample1/outs/filtered_feature_bc_matrix"

# 读取数据
counts <- Read10X(data.dir = data_dir)

# 看一眼矩阵的结构
dim(counts)      # 基因数 × 细胞数
class(counts)    # 应该是 dgCMatrix
counts[1:5, 1:5] # 前5个基因和前5个细胞

注意：如果文件夹里是filtered_feature_bc_matrix.h5这种HDF5格式，那就得用Read10X_h5。我曾在一个项目里把H5文件当普通目录读，折腾了半小时才反应过来。

小提示：Read10X返回的矩阵默认行名为基因ID（ENSEMBL），列名为细胞条形码。如果你想直接拿到基因符号（symbol），可以在读取时加参数gene.column = 2（features.tsv第二列）。不过我习惯后面自己转换，因为有些基因在第二列是空的。

2.2 创建Seurat对象

拿到counts矩阵后，下一步就是把它扔进CreateSeuratObject。这个函数会把矩阵、元数据、技术信息都打包成一个Seurat对象。你想想看，它就像个收纳盒，后面所有的标准化、降维、聚类全在这个盒子里操作。

# 创建Seurat对象
pbmc <- CreateSeuratObject(
  counts = counts,
  project = "PBMC_3k",
  min.cells = 3,    # 至少在3个细胞中表达的基因才保留
  min.features = 200 # 一个细胞至少检测到200个基因
)

# 查看对象
pbmc
# An object of class Seurat 
# 13714 features across 2700 samples within 1 assay 
# Active assay: RNA (13714 features, 0 variable features)

这里min.cells和min.features是我最常用的两个过滤参数。为什么要设它们？说白了就是去掉那些几乎不表达的基因（可能由背景噪音导致）和那些基因数极少的细胞（大概率是死细胞或碎片）。我在处理一个肿瘤样本时，发现min.features设成200还是留了很多低质量细胞，后来改成500才把空泡和双细胞筛干净。具体阈值要根据你的实验来调，别死记硬背。

注意：min.cells和min.features是在创建对象时直接做了一次过滤。如果你不想在这步就过滤，可以设为min.cells = 0, min.features = 0，等后面再用subset手动筛。我个人更建议先保留全部，因为后面我们会用QC指标（线粒体比例、UMI count等）做更细致的筛选。

2.3 数据格式检查——别急着往下跑

对象创建完了，很多人迫不及待就开始标准化。慢着！你得先看看这个对象长什么样，数据有没有读对。我曾经帮同事排查过一个Bug，他跑完NormalizeData结果全是NaN，回来一查——counts矩阵里有一半的细胞列名是重复的。所以，格式检查这步省不得。

我一般会检查以下几项：

检查项	用什么函数	常见问题
细胞数和基因数	`dim(pbmc)`	基因数太少（<500）可能是读取错误
矩阵存储类型	`class(pbmc[["RNA"]]$counts)`	非`dgCMatrix`时后续运算可能变慢
细胞条形码是否唯一	`anyDuplicated(colnames(pbmc))`	重复条形码会导致下游报错
基因名是否有重复	`anyDuplicated(rownames(pbmc))`	重复基因名会影响标准化和差异分析
UMI count分布	`summary(pbmc$nCount_RNA)`	中位数太低（<500）要考虑细胞活性

具体操作看这里：

# 基础维度
dim(pbmc)

# 查看元数据（Seurat自动计算了nCount_RNA和nFeature_RNA）
head(pbmc[[]])

# 检查UMI总数分布
summary(pbmc$nCount_RNA)

# 检查每个细胞表达的基因数
summary(pbmc$nFeature_RNA)

# 看看有没有重复的列名
stopifnot(!anyDuplicated(colnames(pbmc)))

# 看看矩阵是不是稀疏格式
library(Matrix)
is(pbmc[["RNA"]]$counts, "dgCMatrix")

还有一件事：有时候我们从不同来源拿到10X数据，可能features.tsv第三列（类型列）不是"Gene Expression"，或者文件命名不规范。Seurat 4.0以后的Read10X默认会忽略非基因表达的数据，比如抗体捕获（Antibody Capture）或者CRISPR guide。如果你确实需要保留这些，要用Read10X(data.dir, feature.column = 2)之类的参数手动指定。我在一个CITE-seq项目里，抗体标签就存在第四列，当时翻文档翻了好一阵才搞定。

快速自查清单：创建完Seurat对象后，你用下面几行就能完成80%的格式检查。

cat("细胞数:", ncol(pbmc), "\n")
cat("基因数:", nrow(pbmc), "\n")
cat("重复条形码:", anyDuplicated(colnames(pbmc)), "\n")
cat("重复基因名:", anyDuplicated(rownames(pbmc)), "\n")
cat("UMI中位数:", median(pbmc$nCount_RNA), "\n")
cat("基因中位数:", median(pbmc$nFeature_RNA), "\n")

2.4 一张图看懂数据准备流程

下面这张SVG图总结了从10X原始文件到Seurat对象、再到格式检查的完整链路。你可以把它保存下来作为操作路线图。

2.5 常见踩坑与避坑

数据准备阶段我踩过的坑，挑几个典型的说说：

读取H5文件时用错了函数：有的人看到.h5后缀就用Read10X，结果报错cannot open the connection。实际上要换成Read10X_h5。我刚开始也傻傻分不清，直到读了Seurat的vignette才明白。
路径包含中文或空格：Windows用户尤其容易中招。如果路径里有中文，Read10X会读不出来。我用file.choose()交互式选目录来绕过，但更稳妥的办法是养成全英文路径的习惯。
min.cells/min.features设得太严：我有一个样本一共才1000个细胞，把min.cells设成10，结果一半基因都被干掉了。嗯，你要根据数据量调整，别盲目套默认值。
忘记检查细胞条形码重复：有些批次的数据处理时，barcodes.tsv里会有重复行（比如不同泳道合并没去重）。这会导致CreateSeuratObject报错duplicated column names not allowed。解决办法是在读取前用make.unique或自己写个去重脚本。

重点警告：不要在我刚才说的检查步骤跳过去！我见过一个真实案例：有人从10X Cloud拿到的数据，features.tsv里的基因符号是空的（第二列全为NA），结果Read10X读进去行名全是NA，后面所有分析都崩溃。检查时用head(rownames(pbmc))看一眼，能省掉半天调试时间。

2.6 小结

这一章我们走了三个步骤：读取10X数据、创建Seurat对象、检查数据格式。说白了就是“把原料放进锅之前，先看看米有没有发霉”。

我个人习惯把这套检查写成一个小脚本，每次新样本来了跑一遍，比手敲快得多。你可以在自己的项目里也建一个data_ingestion.R文件，把上面的代码揉在一起——第一次配好，后面几乎不用改。