4. 读入空间数据

好，咱们进入实操环节。拿到原始数据后，第一步就是把它读进内存。空间转录组的数据结构比普通单细胞复杂一些，不光是表达矩阵，还有图像、组织坐标、缩放因子。我刚开始接触的时候，就吃过数据路径不对的亏，折腾了半天才发现少读了个文件夹。

这一章我们只讲一件事：如何正确地把10x Visium数据读进来，并创建Seurat或AnnData对象。说白了，就是让数据和代码成功「握手」。

4.1 10x Visium 输出目录结构

先看看10x官方流程跑完后的文件夹长什么样。假设你的样品名叫sample1，典型的输出结构是：

sample1/
├── outs/
│   ├── spatial/
│   │   ├── tissue_positions_list.csv
│   │   ├── scalefactors_json.json
│   │   ├── tissue_hires_image.png
│   │   ├── tissue_lowres_image.png
│   │   └── detected_tissue_image.jpg
│   ├── raw_feature_bc_matrix.h5
│   ├── filtered_feature_bc_matrix.h5
│   └── analysis/
└── sample1_metrics_summary.html

这里面最关键的是两个部分：

• 计数矩阵：filtered_feature_bc_matrix.h5 是官方推荐用的，去掉了空液滴。当然你如果要做质控，也可以从 raw_feature_bc_matrix.h5 开始。
• 空间数据：spatial/ 目录。这是空间转录组特有的，包含组织切片图像、坐标信息和缩放参数。没这个目录，你就没法把表达数据映射回组织图像上。

4.2 读取流程：从文件到 Seurat 对象

在 R 语言里，我们用的是 Seurat 包。从 v4.0 开始，它原生支持空间数据，读起来很方便。看代码：

library(Seurat)
library(patchwork)

# 定义路径
data_dir <- "path/to/sample1/outs/"
spatial_dir <- file.path(data_dir, "spatial")

# 读取计数矩阵 - 关键参数
counts <- Read10X_h5(
  filename = file.path(data_dir, "filtered_feature_bc_matrix.h5"),
  use.names = TRUE  # 保留基因名，不要用 ENSEMBL ID
)

# 创建 Seurat 对象，并添加空间图像
seurat_obj <- CreateSeuratObject(
  counts = counts,
  assay = "Spatial",        # 指定 assay 名为 Spatial
  project = "sample1"
)

# 读取空间图像数据
seurat_obj <- Load10X_Spatial(
  data.dir = data_dir,
  slice = "slice1",
  filename = "filtered_feature_bc_matrix.h5"
)

注意：Read10X_h5 和 Load10X_Spatial 是两个不同的函数。前者只读表达矩阵，后者会把 spatial/ 目录下的图像和坐标一起读进来。如果你用前者，还得手动加图像，比较麻烦。我建议直接用 Load10X_Spatial，一步到位。

4.3 读取流程：从文件到 AnnData 对象

Python 这边，我们用 Scanpy 生态。不过 Scanpy 本身没有空间数据模块，需要搭配 squidpy 或 scimap。我个人喜欢用 squidpy，它和 Scanpy 配合得很默契。

import scanpy as sc
import squidpy as sq
import numpy as np

# 定义路径
data_dir = "path/to/sample1/outs/"
adata = sq.read.visium(
    path = data_dir,
    library_id = "sample1",
    load_images = True,      # 是否加载组织图像
    source_image_path = None # 如果要用高分辨率图，填路径
)

# 看一下数据
print(adata)
print(adata.obsm['spatial'])
print(adata.uns['spatial'])

sq.read.visium 这个函数会帮你完成三件事：读取 h5 文件、解析 tissue_positions_list.csv、加载图像缩放参数。返回的 annData 对象里：

• adata.X：表达矩阵。
• adata.obsm['spatial']：每个 spot 的 tissue 坐标（行坐标、列坐标）。
• adata.uns['spatial']：图像路径和缩放因子。

如果你用 Scanpy 自带的 sc.read_visium，也是类似的。不过我记得 sc.read_visium 在早期版本里会少读一些元数据，所以我还是推荐 sq.read.visium。

4.4 核心数据结构对比

我把 Seurat 和 AnnData 在空间数据上的对应关系整理了一下。理解了这张表，你在两个生态之间切换就不会迷路。

4.5 常见问题与避坑

读了这么多次空间数据，我总结出几个高频报错点：

1. h5 文件版本不兼容：10x 偶尔会更新 h5 的内部格式。如果报错说「不支持的 h5 版本」，建议把 Seurat 或 Scanpy 升级到最新版。
2. tissue_positions_list.csv 列名不一致：不同版本的 spaceranger 输出的列名可能不一样。有的版本用 in_tissue，有的用 tissue。读进来之后，记得先检查一下列名。
3. 图像路径混乱：如果你把数据复制到别的机器上，图像文件路径可能变成绝对路径。这时候需要手动修正路径，或者重新读一遍。
4. 内存不够：高分辨率组织图像（tissue_hires_image.png）可能很大。如果你只做分析不需要看图，可以在读数据时把 load_images = False，节省内存。

4.6 读取后的验证

数据读进来了，怎么知道它对不对？我一般做三个快速检查：

• 检查 spot 数量：filtered 数据一般有 3000-5000 个 spot，具体取决于组织大小。如果读出来只有几十个，大概率是读错了 reference。
• 检查基因数量：Visium 捕获大约 20000 个基因。如果少很多，可能用了 wrong 的 h5 文件。
• 检查坐标范围：空间坐标应该在几百到几千的范围内，和图像像素对应。

# Seurat 快速检查
ncol(seurat_obj)  # spot 数
nrow(seurat_obj)  # 基因数
head(GetTissueCoordinates(seurat_obj))

# AnnData 快速检查
adata.n_obs   # spot 数
adata.n_vars  # 基因数
adata.obsm['spatial'][:5, :]

如果这三个都没问题，恭喜你，数据已经成功读进来了。接下来就可以开始做批次质控了。

4.7 多切片读取准备

如果你的项目涉及多个切片（比如同一个病人的不同区域，或者不同病人的切片），建议从最开始就把数据组织规范。我一般是建一个列表，每个元素是一个样品：

# Seurat 多切片
sample_dirs <- c("sample1/outs/", "sample2/outs/", "sample3/outs/")
seurat_list <- list()

for (i in seq_along(sample_dirs)) {
  seurat_list[[i]] <- Load10X_Spatial(
    data.dir = sample_dirs[i],
    slice = paste0("slice", i)
  )
}

# 合并
merged_obj <- merge(seurat_list[[1]], y = seurat_list[-1])

# AnnData 多切片 - 建议用 squidpy
adatas = {}
for sid in ["sample1", "sample2", "sample3"]:
    adatas[sid] = sq.read.visium(f"data/{sid}/outs/", library_id=sid)
    adatas[sid].obs['batch'] = sid

# 后续可以用 sc.concat 合并
adata_all = sc.concat(adatas, label="batch", index_unique="-")

注意：合并之前，先确保每个切片的基因名一致。10x 使用不同的参考基因组版本时，基因名可能不同（比如 ENSG000... 和 GENE_NAME）。不一致的话，合并后会损失很多基因。