📚 课程目录第4章空间转录组数据准备

空间转录组数据准备：10x Visium数据导入与空间坐标

好，开始咱们这一章。上一章我们把 CellTrek 的环境搭好了，但巧妇难为无米之炊，对吧？这一章咱们就把“米”——也就是空间转录组数据——给准备好。

我直接用 10x Visium 的数据来做演示。为什么选它？因为它是目前最主流的平台，数据格式规范，文档也全。你先跟它混熟了，其他平台的数据处理思路其实大同小异。

10x Visium 数据长什么样？

你从 10x 官网下载或者自己跑完 Space Ranger 后，会拿到一个目录。里面最重要的就是下面这三个文件。我把它们整理成一个表格，你一眼就能看明白。

文件/目录名	说明	我的使用习惯
filtered_feature_bc_matrix.h5	过滤后的基因表达矩阵 (HDF5格式)	我一般直接用这个，省去读三个文件的麻烦
tissue_positions_list.csv	每个捕获点的空间坐标	这个文件必须仔细检查，坐标格式有时会微调
scalefactors_json.json	图像缩放参数、组织半径等	CellTrek 不直接用它，但我用来做坐标验证

核心思路：我们要做的，就是把上面这堆文件里的信息，合并成一个 Seurat 对象。这个对象里既有表达数据，又有空间坐标。

第一步：把基因表达数据读进来

我个人习惯用 Read10X_h5() 这个函数。它一步到位，直接把 h5 文件读成矩阵。比你用 Read10X() 读三个 mtx 文件要清爽得多。

# 加载Seurat包
library(Seurat)

# 定义数据路径
data_dir <- "path/to/your/visium_data/"

# 读取过滤后的h5文件
counts <- Read10X_h5(paste0(data_dir, "filtered_feature_bc_matrix.h5"))

# 看一眼矩阵长啥样
dim(counts)
# 一般会输出：特征数 x 捕获点数

这里有个小细节需要注意。你读进来的 counts 默认是一个 dgCMatrix 格式的稀疏矩阵。为什么要用稀疏矩阵？因为空间转录组数据里，绝大多数基因在单个捕获点里是检测不到的，0 值特别多。稀疏矩阵能省下 80% 以上的内存。我刚开始做的时候没注意这点，直接用普通矩阵，结果 8G 内存的机器直接崩了。

第二步：读入并解析空间坐标

空间坐标都藏在 tissue_positions_list.csv 这个文件里。嗯，这里要注意，它的格式可能让你有点懵。咱们先把它读进来看看。

# 读入空间坐标文件
spatial_coords <- read.csv(paste0(data_dir, "spatial/tissue_positions_list.csv"),
                           header = FALSE,
                           row.names = 1)

# 给列名起个好记的名字
colnames(spatial_coords) <- c("in_tissue", "array_row", "array_col", 
                               "pxl_row_in_fullres", "pxl_col_in_fullres")

# 只看前几行
head(spatial_coords)

你拿到手的表格应该是这样的。

行名 (barcode)	in_tissue	array_row	array_col	pxl_row_in_fullres	pxl_col_in_fullres
AAACAAGTATCTCCCA-1	1	0	16	512	1024
AAACACCAATAACTGC-1	1	0	24	768	512
AAACCCGAACGAAATC-1	0	1	4	...	...

这里的关键列是 in_tissue。它告诉你这个捕获点是不是在组织覆盖的区域上。1 表示是，0 表示不是。我们在后续分析中，一般只保留 in_tissue == 1 的点。那些在玻片背景上的点，表达量极低，留着只会增加噪声。

实用技巧：你在筛选的时候，用 row.names(spatial_coords[spatial_coords$in_tissue == 1, ]) 就能拿到所有有效点的 barcode 列表。然后用这个列表去过滤你的表达矩阵。

第三步：合并数据，构建 Seurat 对象

数据有了，坐标也有了。接下来就是把它们捏合在一起。我建议用 CreateSeuratObject() 先创建基础对象，然后用 AddMetaData() 把坐标信息添加上去。

# 先创建一个基础的Seurat对象
# 只保留组织上的点
tissue_barcodes <- rownames(spatial_coords[spatial_coords$in_tissue == 1, ])
counts_filtered <- counts[, tissue_barcodes]

brain_obj <- CreateSeuratObject(counts = counts_filtered,
                                 project = "Visium_Brain",
                                 assay = "Spatial")

# 把坐标信息加到meta.data里
brain_obj <- AddMetaData(brain_obj, 
                          metadata = spatial_coords[tissue_barcodes, ])

# 确认一下
head(brain_obj[[]])

好了，现在 brain_obj 里就包含了表达数据和坐标信息。不过，CellTrek 要求坐标是标准化的，单位是微米。而目前我们取的坐标是像素值。所以还需要一步转换。

至于转换公式，10x 官方文档里提到，1 个 spot 的直径对应的像素数是 55 个像素。如果你需要把像素坐标转成微米坐标，可以用下面的方法算个近似值。不过在实际做 CellTrek 的时候，我发现直接用像素坐标，然后让 CellTrek 自己根据点间距离来推断组织空间结构，效果往往更好。这点我们下一章再细聊。

避坑指南：我曾经在这个地方踩过坑。当时我按照网上一些教程，把坐标做了归一化，又做了中心化，结果 CellTrek 跑出来的结果全乱套了。后来我发现，CellTrek 在内部会自己做坐标的标准化处理。你传给它原始像素坐标就行，千万别自己画蛇添足。

第四步：快速验证数据完整性

数据构建完成后，我强烈建议你做两个最基础的检查。别急着往下跑，花两分钟看一眼，能省下后面几个小时的调试时间。

检查一：细胞数和基因数是否合理

# 查看对象的基本信息
brain_obj
# 应该能看到类似：14240 features across 3003 samples
# 如果只有几十个细胞，那可能筛选的时候把组织点滤掉了

检查二：坐标范围是否合理

# 看看坐标的分布范围
summary(brain_obj$pxl_row_in_fullres)
summary(brain_obj$pxl_col_in_fullres)

# 一般 Visium 的捕获区域是 6.5mm x 6.5mm，
# 对应的像素坐标大概在 0~2000 之间。
# 如果出现负数或超大值，那肯定是读错了。

经验之谈：我遇到过好几次这样的情况——数据看起来读进来了，坐标也对，但跑 CellTrek 时老是报错。后来发现是 barcode 名字的格式不统一。表达矩阵的列名是 AAACAAGTATCTCCCA-1，而坐标文件的行名是 AAACAAGTATCTCCCA_1。一个短横一个下划线，坑了我两个小时。建议你读进来后，用 gsub("-", "_", ...) 统一一下格式。

好了，数据准备就到这儿。现在你的 R 环境里已经有了一个结构完整的 Seurat 对象，而且空间坐标信息也附在 meta.data 里了。这个对象就是我们后续所有分析的基础。