4. 空间转录组数据格式：h5ad、STARmap、10x Visium 解析、读取与基础操作（使用 scanpy）

好，咱们开始干正事。上一章聊了聊空间转录组的基本概念，你可能已经跃跃欲试想上手操作了。但在此之前，有件事必须搞清楚——数据格式。

在我刚入行那几年，最头疼的就是数据格式不统一。各家平台各玩各的，h5ad、HDF5、JSON、CSV、TIFF……每次换个平台做分析，光读数据就得花半天。后来社区逐渐收敛，h5ad 成了事实标准，但 STARmap 和 10x Visium 的数据组织方式还是大不相同。这章我就带着你把这些格式的门道理清楚，咱们动手读一读，心里才踏实。

4.1 格式概览：它们到底存了什么？

说白了，空间转录组数据就是“一张地图 + 每个坐标上的基因表达”。不同平台的区别在于：地图是怎么画的，表达矩阵怎么放的。

我个人习惯先看数据结构的抽象层。下面这张图可以帮你快速建立认知框架：

你看，虽然源头格式五花八门，但最终都会被 scanpy 处理成统一的 AnnData 对象。这正是 scanpy 的价值所在——它帮你屏蔽了底层格式差异。

4.2 h5ad 格式：scanpy 的母语

h5ad 是 scanpy 的原生格式，基于 HDF5 构建。你可以把它理解成一个超级容器：基因表达矩阵、细胞元数据、基因元数据、空间坐标、甚至图像，全都能塞进去。

我最推荐的做法是：不管数据从哪来，最终都存成 h5ad。这样下次再用就不用折腾了。

4.2.1 读取 h5ad

读取 h5ad 文件，一行代码就够了：

import scanpy as sc

adata = sc.read_h5ad('my_data.h5ad')
print(adata)

输出通常是这样的：

AnnData object with n_obs × n_vars = 5000 × 20000
    obs: 'cell_type', 'batch', 'n_genes_by_counts'
    var: 'gene_symbols', 'highly_variable'
    uns: 'spatial'
    obsm: 'spatial'

看到 obsm: 'spatial' 了吧？这就是空间坐标的存储位置。在 AnnData 里，obsm 专门放“观测 x 特征”的辅助矩阵，空间坐标就放在这里。

4.2.2 探究内部结构

你拿到一个 h5ad 文件，第一件事是什么？我习惯先看看里面到底有啥：

# 看看 obs 里有哪些列
print(adata.obs.columns.tolist())

# 看看 var 里有哪些列
print(adata.var.columns.tolist())

# 看看 obsm 里有哪些键
print(list(adata.obsm.keys()))

# 看看 uns 里有什么
print(list(adata.uns.keys()))

4.2.3 构建自己的 h5ad

如果你有一堆数据想整理成 h5ad，也很简单：

import numpy as np
import pandas as pd

# 假设你准备好了这些
X = np.random.rand(100, 5000)  # 100个细胞，5000个基因
cell_metadata = pd.DataFrame({
    'cell_id': range(100),
    'region': ['cortex'] * 100
})
gene_metadata = pd.DataFrame({
    'gene_name': [f'Gene_{i}' for i in range(5000)]
})
coordinates = np.random.rand(100, 2)  # 空间坐标

# 构建 AnnData
from anndata import AnnData
adata = AnnData(X=X, obs=cell_metadata, var=gene_metadata)
adata.obsm['spatial'] = coordinates

# 保存
adata.write('my_spatial_data.h5ad')

嗯，这里要注意：obs 的索引要跟 X 的行一一对应；var 的索引要跟 X 的列一一对应。

4.3 STARmap 格式解析

STARmap 的数据格式有点特殊。它不像 10x 那样把所有数据打包成几个标准文件，而是以多个文件形式存在。我在项目中遇到过好几次，搞混文件名就全废了。

4.3.1 典型文件结构

一个标准的 STARmap 数据集通常包含：

文件	内容
cell_barcode_count.txt 或 *.csv	细胞 x 基因 表达矩阵
cell_metadata.csv	细胞 ID、类型、区域等元信息
cell_position.csv	空间坐标信息（x, y, z）
images/ 文件夹	原始图像栈（可选）

4.3.2 读取流程

我的标准做法是三步走：先读矩阵，再读元数据，最后整合到一起：

import pandas as pd
import scanpy as sc

# 1. 读表达矩阵
counts_df = pd.read_csv('cell_barcode_count.csv', index_col=0)
# 确保是细胞名在行、基因名在列

# 2. 读空间坐标
pos_df = pd.read_csv('cell_position.csv', index_col=0)
# 通常有 'x', 'y' 两列

# 3. 读元数据
meta_df = pd.read_csv('cell_metadata.csv', index_col=0)

# 4. 整合成 AnnData
adata = sc.AnnData(counts_df)
adata.obs = meta_df
adata.obsm['spatial'] = pos_df[['x', 'y']].values

print(adata)

4.3.3 对齐检查技巧

我一般会加一个验证：

# 检查索引是否一致
print("矩阵细胞数量:", counts_df.shape[0])
print("坐标细胞数量:", pos_df.shape[0])
print("元数据细胞数量:", meta_df.shape[0])

# 检查索引是否完全相同
print("矩阵 vs 坐标:", counts_df.index.equals(pos_df.index))
print("矩阵 vs 元数据:", counts_df.index.equals(meta_df.index))

# 如果不一致，需要按交集对齐
common_cells = counts_df.index.intersection(pos_df.index)
common_cells = common_cells.intersection(meta_df.index)

counts_df = counts_df.loc[common_cells]
pos_df = pos_df.loc[common_cells]
meta_df = meta_df.loc[common_cells]

4.4 10x Visium 格式解析

10x Visium 是目前最常用的空间转录组平台，它的数据格式相对标准化。你下载完数据后，通常看到一个文件夹结构：

spaceranger_output/
├── filtered_feature_bc_matrix.h5
├── filtered_feature_bc_matrix/
│   ├── barcodes.tsv.gz
│   ├── features.tsv.gz
│   └── matrix.mtx.gz
└── spatial/
    ├── tissue_positions_list.csv
    ├── scalefactors_json.json
    ├── tissue_hires_image.png
    └── tissue_lowres_image.png

这里有两个关键文件：filtered_feature_bc_matrix.h5（HDF5 格式的表达矩阵）和 spatial/tissue_positions_list.csv（spot 坐标）。

4.4.1 标准读取方法

scanpy 自带的 read_visium 函数就是为这个准备的：

adata = sc.read_visium(
    'spaceranger_output/',   # 10x 输出目录的路径
    count_file='filtered_feature_bc_matrix.h5',  # 可以用 h5 或 mtx
    library_id='sample_001',  # 给样本起个名字
    load_images=True          # 是否加载组织图像
)

print(adata)

输出跟我之前展示的类似，但多了几个特殊字段：

AnnData object with n_obs × n_vars = 4000 × 25000
    obs: 'in_tissue', 'array_row', 'array_col'
    var: 'gene_symbols', 'feature_types'
    uns: 'spatial'
    obsm: 'spatial', 'spatial_original'

4.4.2 手工读取（不依赖 read_visium）

如果你装 scanpy 版本比较旧，或者想深入理解数据组织，也可以手动读：

import anndata
import pandas as pd
import scipy.io

# 读 barcodes
barcodes = pd.read_csv(
    'filtered_feature_bc_matrix/barcodes.tsv.gz',
    header=None, sep='\t'
).iloc[:, 0].tolist()

# 读 features
features = pd.read_csv(
    'filtered_feature_bc_matrix/features.tsv.gz',
    header=None, sep='\t'
).iloc[:, 1].tolist()  # 第二列是基因名

# 读矩阵
mtx = scipy.io.mmread(
    'filtered_feature_bc_matrix/matrix.mtx.gz'
).tocsr()

# 读空间坐标
positions = pd.read_csv(
    'spatial/tissue_positions_list.csv',
    header=None,
    names=['barcode', 'in_tissue', 'array_row', 'array_col', 'pxl_col_in_fullres', 'pxl_row_in_fullres']
)

# 构建 AnnData
adata = anndata.AnnData(X=mtx.T, obs=pd.DataFrame(index=barcodes), var=pd.DataFrame(index=features))

# 对齐坐标
pos_idx = positions.set_index('barcode')
adata.obs['in_tissue'] = pos_idx.loc[barcodes, 'in_tissue'].values
adata.obsm['spatial'] = pos_idx.loc[barcodes, ['pxl_col_in_fullres', 'pxl_row_in_fullres']].values

print(adata)

手动走一遍，你就能彻底理解 scanpy 的 read_visium 在背后做了什么。

4.5 用 scanpy 做点基础操作

数据读进来了，总得干点正事。我挑三个最基础、最常用的操作演示给你看。

4.5.1 查看与切片

拿到 AnnData 对象后，你首先要学会怎么看它、怎么切它：

# 查看基本统计
print(adata.obs['n_genes_by_counts'].describe())

# 选择 in_tissue 的 spot
adata_tissue = adata[adata.obs['in_tissue'] == 1].copy()

# 选择某些基因
adata_subset = adata[:, ['GeneA', 'GeneB', 'GeneC']]

# 查看某个基因的表达分布
print(adata[:, 'GeneA'].X.toarray().flatten()[:10])

切完之后记得 .copy()。AnnData 默认是引用，改切片可能会改原数据。这个坑我踩过不止一次。

4.5.2 过滤低质量 spot/细胞

扫一眼数据质量，过滤掉那些表现太差的：

sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)

# 手动过滤：去除高线粒体比例
adata.var['mt'] = adata.var_names.str.startswith('MT-')
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False)
adata = adata[adata.obs['pct_counts_mt'] < 20, :].copy()

4.5.3 归一化与对数变换

归一化是所有定量分析的前提：

# 总 UMIs 归一化
sc.pp.normalize_total(adata, target_sum=1e4)

# 对数变换
sc.pp.log1p(adata)

# 高变基因选取（加速后续分析）
sc.pp.highly_variable_genes(adata, n_top_genes=2000)
adata = adata[:, adata.var['highly_variable']].copy()

我个人习惯把归一化和对数变换连起来做。很多新手会忘记对数化，结果 PCA 一跑全是一个方向——表达量差异太大。

4.6 三种格式的快速对比

我把核心区别整理成了一张表，方便你对照：

维度	h5ad	STARmap	10x Visium
格式类型	单文件（HDF5）	多文件（CSV/txt）	多文件（HDF5 + 图像 + CSV）
读取方式	sc.read_h5ad()	手动组装	sc.read_visium()
空间坐标位置	obsm['spatial']	单独 CSV 文件	spatial/tissue_positions_list.csv
图像嵌入	可选（h5ad 内）	单独图像文件夹	默认嵌入 uns 中
社区主推	是（scanpy 原生）	否	部分（有官方支持）
推荐做法	最终统一转 h5ad	读进来后存 h5ad	用 read_visium 直接读