📚 课程目录第3章高效数据存储格式

3、高效数据存储格式：HDF5、Zarr、TileDB在空间转录组中的应用

空间转录组数据有多大？我给你们算笔账。

一张10x Visium切片，大概5000个spots，每个spot测到几千个基因。原始数据出来，少说几个GB。要是换成Slide-seq或者MERFISH这种高分辨率的，几十GB甚至上百GB都很正常。你想想看，一整个项目下来，几十张切片，数据量轻松上TB。传统的小文件、CSV格式？嗯，基本跑不动。

所以，存储格式的选择，直接决定了你的分析流程是「等几分钟」还是「等到天荒地老」。我这些年折腾下来，真正能在生产环境扛住空间转录组数据压力的，就三个：HDF5、Zarr、TileDB。今天我把它们掰开揉碎了讲清楚。

3.1 HDF5：经典方案，但要注意坑

HDF5，学术界的「老前辈」了。单细胞领域的主力Scanpy，它的底层存储就是HDF5格式（.h5ad文件）。我刚开始做空间转录组时，也直接从它入手，图个省事。

HDF5的本质，是一个文件系统的抽象。里面可以创建「组」（类似文件夹），组里面放「数据集」（类似文件）。你可以在同一个.h5文件中，把表达矩阵、空间坐标、元数据全部塞进去，结构清晰。

典型结构：一个.h5ad文件内部看起来像这样——

/X            # 主表达矩阵 (shape: n_obs x n_vars)
/obs          # 细胞/spot 元数据表
/var          # 基因元数据表
/obsm         # 降维结果或空间坐标
/uns          # 非结构化数据（聚类marker等）

当初我接手的第一个Visium项目，数据量大概200GB，全用HDF5存。跑Leiden聚类时，发现加载一个完整矩阵就要等5分钟。为什么会这样？因为HDF5的读写是「序列化」的，你读大矩阵时，它要先把整个数据块从磁盘搬到内存，再解压，再切片。

警告：HDF5不是不能做并行，但需要MPI支持，配置起来相当麻烦。我在集群上折腾过两回，每次都被MPI版本兼容问题搞得头大。如果你的计算节点是多核并行读取同一个HDF5文件，强烈建议先做压力测试，否则容易遇到锁竞争导致的性能崩塌。

HDF5的优势：工具链极其成熟。R语言的rhdf5、Python的h5py、MATLAB内置支持，几乎零门槛上手。加上强大的数据压缩能力（比如搭配gzip或LZF），磁盘占用能省不少。

HDF5的短板：单文件设计，并行写入就是噩梦；而且文件损坏后恢复困难。我记得有一次服务器意外断电，一个写了12小时的HDF5文件直接打不开了，那叫一个心痛。从那以后，我但凡跑长时间处理任务，都会在代码里加检查点备份。

3.2 Zarr：云原生的新贵

Zarr的设计哲学，跟HDF5正好相反。HDF5追求「把一切塞进一个文件」，Zarr追求「把数据切碎，每个碎片独立存储」。说白了，Zarr把一个大数组切成了许多小块（chunk），每个块单独压缩、单独存储。可以存在本地文件夹，也可以直接存在S3或GCS上。

个人习惯：我在做大规模MERFISH数据分析时，全部改用Zarr。原因很简单——我可以在云端开100个worker，同时去读同一个Zarr数据集的不同切片，完全不会互相干扰。这在HDF5里想都不敢想。

Zarr的块存储带来了几个颠覆性的好处：

并行读取无敌：每个chunk是一个独立文件，多线程/多进程读取天然安全。
惰性加载：只需要读取触及到的chunk，其余根本不会碰到磁盘。体现在代码层面就是，你只取一小块数据，读的速度极快。
云原生：存在S3上的Zarr数据集，可以直接用zarr.open()远程读取，不需要先下载。我有个项目，数据放在阿里云OSS上，分析脚本跑在本地服务器，数据流全程不走本地磁盘。

看一个简单的例子，用Zarr创建一个空间表达矩阵并写数据：

import zarr
import numpy as np

# 创建一个10000个spot x 2000个基因的压缩数组
store = zarr.open('spatial_data.zarr', mode='w',
                  shape=(10000, 2000),
                  chunks=(1000, 200),
                  dtype='float32',
                  compressor=zarr.Blosc(cname='zstd', clevel=3, shuffle=2))

# 写入模拟数据
data = np.random.rand(10000, 2000).astype('float32')
store[:] = data

# 读取某一个spot的基因表达
spot_426 = store[426, :]

  补充细节：上面代码中，chunks=(1000, 200) 意味着每个块是1000行×200列。chunk的大小直接影响性能。太小了导致文件碎片过多，太大了又失去惰性加载优势。我的经验是，chunk大小控制在1~10MB之间比较合适，既能压缩好，又不会过多碎片。

Zarr的局限：它没有原生的事务写入支持。如果你两个进程同时写同一个Zarr数组的不同区域，理论上可以，但实际操作中得小心别覆盖元数据。我也踩过这个坑，后来改用TileDB解决并发写入问题。

3.3 TileDB：生产级的多维数组引擎

TileDB在单细胞社区里可能没前两者那么普及，但它在生信生产环境中越来越受重视。本质里，TileDB是一个多维数组的存储引擎，支持密集数组和稀疏数组两种模式。

空间转录组数据天然适合稀疏数组——因为大部分基因在单个spot里是0。TileDB的稀疏数组存储只会记录非零值的位置和数值，省空间离谱。我做过一个测试，同样一份10x Visium的原始UMI计数矩阵，CSV格式是3.2GB，压缩成TileDB稀疏数组后只有120MB。

TileDB的杀手特性：

ACID事务：支持原子写，不怕写到一半断电。我从前做数据管道时最头疼的就是「写到一半挂了，整个文件废了」，用TileDB之后彻底没这个烦恼。
过滤与聚合下推：可以在读取时就做条件过滤，比如「只返回表达量大于5的基因」，这部分计算在存储引擎层面完成，不用把整个数组先拉到内存。
多语言绑定：C++、Python、R、Java、甚至Julia都有官方绑定。我在一个C++编写的空间分析工具里直接嵌入了TileDB作为后端，性能比用Python中转快了一个数量级。

一个简单的Python例子，用TileDB存储空间坐标和表达：

import tiledb
import numpy as np

# 创建稀疏数组描述
ctx = tiledb.Ctx()
schema = tiledb.ArraySchema(
    domain=tiledb.Domain(
        tiledb.Dim(name='spot_id', domain=(0, 9999), tile=100),
        tiledb.Dim(name='gene_id', domain=(0, 1999), tile=50)
    ),
    attrs=[tiledb.Attr(name='expression', dtype=np.float32)],
    sparse=True,
    ctx=ctx
)

# 写入数据
with tiledb.SparseArray('spatial_tiledb', mode='w', schema=schema) as arr:
    # 只写入非零值
    spots = np.array([0, 0, 1, 1, 2], dtype=np.uint64)
    genes = np.array([10, 50, 30, 100, 200], dtype=np.uint64)
    vals = np.array([8.0, 3.5, 12.0, 6.2, 9.1], dtype=np.float32)
    arr[spots, genes] = vals

# 读取spot 0的所有非零表达
with tiledb.SparseArray('spatial_tiledb', mode='r') as arr:
    result = arr.query(cond='spot_id == 0')[:]
    print(result)

TileDB要注意的地方：它引入了自己的概念，学习曲线比Zarr稍陡。另外，TileDB的数组是静态schema设计（维度、属性在创建时固定），如果你需要频繁修改整体结构，得慎重。我的建议是：一旦设计好schema，尽量不要改，避免重写大量数据。

3.4 横向对比：怎么选？

我倾向于这样判断：

维度	HDF5	Zarr	TileDB
文件结构	单一文件	目录/云存储碎片	目录/云存储碎片
并行读取	弱（需MPI）	强（天然分块）	强（原生并行）
并行写入	极弱（单一文件锁）	中等（块级独立）	强（ACID事务）
云原生程度	低	高	高
压缩与效率	好（gzip/LZF）	极好（Blosc）	极好（多种可选）
学习成本	低	低	中
适用范围	小团队/单机工具	云分布式/大规模	生产管道/数据仓库

  避坑指南：我曾经在一个项目中，场景是「多个工作节点同时写一份中间结果」，当时选型了Zarr，结果频繁出现metadata冲突。最后换成TileDB才解决。反过来，如果是做「一次写入、多次只读」的分析，Zarr的便利性和性能优势完全足够，没必要上TileDB增加复杂度。

总结一下，不做绝对的好坏评判，但有一条铁律：根据数据规模、计算场景和团队技术栈来选，不要因为「大家都在用」就盲目跟风。我经手的每个项目，选型前都会做一个小规模的压力测试，模拟真实读写负载，这比任何理论分析都靠谱。