📚 课程目录第3章单细胞RNA-seq数据准备

单细胞RNA-seq数据准备：格式要求、Seurat对象创建

数据准备是CellTrek流程中最基础、也最容易翻车的一步。说白了，你后面所有漂亮的可视化结果都建立在数据格式正确的基础上。我见过不少同学在这步上栽跟头——读进去的数据是乱码，或者Seurat对象死活建不起来，浪费大量时间排查。

今天我们就来把这块彻底捋清楚。看完这章，你必须能独立完成从原始表达矩阵到标准Seurat对象的构建，并且知道每种格式的坑在哪里。

一、输入格式：你拿到的数据长什么样？

单细胞测序的下游数据五花八门，但最主流的还是10X Genomics平台输出的格式。我个人习惯先把原始数据过一遍，看看文件夹结构，再做决定。

1.10X Genomics官方格式（mtx + barcodes + features）

这是CellTrek最推荐的形式。你会得到三个文件：

matrix.mtx：Market Exchange Format，稀疏矩阵，存的是UMI计数
barcodes.tsv：每个细胞条的ID
features.tsv：基因名（旧版叫genes.tsv）

小技巧：我排查问题时第一件事就是检查这三个文件的列数是否匹配。尤其是features.tsv，现在Seurat要求至少两列（gene_id和gene_symbol），缺一列就会报错。

2.H5格式（Cell Ranger v3+）

一个HDF5文件搞定所有，基因信息、条形码、矩阵全包在里面。好处是文件少，但你需要确认H5是v3格式还是更老的v2格式。之前有学员拿了个v2的H5文件，Read10X_h5()读了却报维度错误——后来发现Seurat 4.x对v2支持不太友好，得手动改版本。

# 推荐使用Read10X_h5()读取v3格式
library(Seurat)
counts <- Read10X_h5("filtered_feature_bc_matrix.h5")
dim(counts)

3.CSV/TSV表格

有时候上游分析给了你一个表达矩阵CSV，行是基因，列是细胞。这种格式最容易出问题——分隔符不统一、行列号缺失、基因名有重复。我举个例子：

# 读取CSV为稀疏矩阵（注意行为基因，列为细胞）
data <- read.csv("expr_matrix.csv", row.names = 1, check.names = FALSE)
# 检查行列名
head(rownames(data))  # 应该都是基因
head(colnames(data))  # 应该都是细胞ID
# 转为稀疏矩阵（节省内存）
library(Matrix)
sparse_mat <- as(as.matrix(data), "dgCMatrix")

注意：CSV文件里的基因名如果包含特殊字符（比如“-”、“.”），read.csv默认会帮你改成合法R名称。你应该用check.names = FALSE避免被篡改。我吃过这个亏，跑CellTrek时基因名对不上，排查了整整一上午。

4.HDF5 / AnnData（Python生态）

如果你用的是Scanpy处理过的.h5ad文件，可以用SeuratDisk包转换。但我更推荐直接读原始10X产物，因为转换过程中可能会丢元数据。

# 将h5ad转成Seurat（需要SeuratDisk）
library(SeuratDisk)
Convert("data.h5ad", dest = "h5seurat", overwrite = TRUE)
seurat_obj <- LoadH5Seurat("data.h5seurat")

二、创建Seurat对象：把数据装进容器

数据拿到手，下一步就是用CreateSeuratObject()把它包装成一个标准的Seurat对象。这一步非常直白，但细节决定成败。

核心命令

# 从10X的mtx文件夹创建
counts <- Read10X(data.dir = "path/to/filtered_feature_bc_matrix/")
seurat_obj <- CreateSeuratObject(counts = counts, 
                                  project = "CellTrek_Demo", 
                                  min.cells = 3, 
                                  min.features = 200)

上面min.cells = 3表示某个基因至少在3个细胞中被检测到才保留；min.features = 200则要求每个细胞至少有200个基因。这两个参数是经验值，但你在实际项目中可能需要调整。

个人经验：如果不做预过滤就直接创建对象，后面跑QC时会看到一堆全是零的垃圾细胞，徒增计算量。我一般都会设置min.features = 200以上，除非数据质量特别差。

添加样本元数据

很多时候你有多个样本，需要为每个细胞打上样本标签。在创建对象之后，用AddMetaData()加进去：

# 假设barcodes向量对应每个细胞的样本ID
sample_info <- data.frame(barcode = colnames(seurat_obj), 
                           sample = rep("Sample1", ncol(seurat_obj)))
rownames(sample_info) <- sample_info$barcode
seurat_obj <- AddMetaData(seurat_obj, metadata = sample_info)

三、质量检查与清洗：扔掉坏细胞

创建好对象后，千万别急着跑下游分析。先看看数据质量怎么样。我习惯用三个指标过滤：

指标	含义	典型过滤范围	注意事项
nFeature_RNA	检测到的基因数量	200 ~ 3000	过高可能是双细胞/空液滴
nCount_RNA	总UMI数量	500 ~ 10000	与nFeature正相关，离群值需检查
percent.mt	线粒体基因比例	< 20% (人类) / < 10% (小鼠)	过高说明细胞正在死亡

计算线粒体比例需要一点小技巧：

# 计算线粒体基因百分比
seurat_obj[["percent.mt"]] <- PercentageFeatureSet(seurat_obj, pattern = "^MT-")

然后画个提琴图看看分布：

VlnPlot(seurat_obj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

过滤操作：

seurat_obj <- subset(seurat_obj, subset = nFeature_RNA > 200 & nFeature_RNA < 3000 & percent.mt < 20)

注意：过滤后最好再画一次图，看看是不是仍有离群值。

避坑指南：我曾经在一个免疫组数据里发现线粒体比例高达60%，以为是细胞死亡，后来才发现样本处理过程中用的组织消化酶导致了RNA降解。此时不应简单过滤，而要考虑是否使用Hashtag或细胞标签去除死细胞。数据异常先想原因，不要无脑过滤。

四、标准化：让细胞之间可比较

过滤完，下一步是标准化。CellTrek要求输入归一化后的数据，这里用最常见的LogNormalize：

seurat_obj <- NormalizeData(seurat_obj, normalization.method = "LogNormalize", scale.factor = 10000)

这个步骤把每个细胞的UMI总数归一化到10000，然后取log(1+x)。经过这一步，基因表达量就可以跨细胞比较了。

最终整理：一张表总结创建流程

步骤	函数	输出	备注
1. 读取数据	Read10X / read.csv	稀疏矩阵	注意列顺序
2. 构建对象	CreateSeuratObject	Seurat对象	设好min.cells/min.features
3. 计算QC	PercentageFeatureSet	meta.data新增列	线粒体基因用pattern
4. 过滤	subset	过滤后的Seurat对象	根据nFeature与MT%
5. 标准化	NormalizeData	标准化后的Seurat对象	默认LogNormalize

到这步，你的Seurat对象就达到了CellTrek的数据准入标准。记住：输入格式、对象创建、QC过滤、标准化，这四个环节缺一不可。我建议你在每次创建完对象后，都跑一个str(seurat_obj)确认结构完整，特别是assays$RNA@counts和meta.data列名。

实操中如果遇到问题，优先检查：文件路径是否有空格？barcodes.tsv有没有空行？CSV分隔符是逗号还是制表符？这些看似低级的问题，其实占了我80%的排错时间。