4、稀疏矩阵表示：用 scipy.sparse 把内存省下来

空间单细胞数据，说白了就是个天然稀疏的东西。一个 spot 里可能只有几百到几千个基因表达量非零，剩下的几万个全是零。如果你老老实实存成稠密矩阵，那内存简直是在烧钱。我刚入行那会儿，手头一个 10× 空间切片，基因数两万，spot 数五千，用 np.float32 一算——400 MB，还没做 QC 就快把笔记本撑爆了。后来換成 CSR 格式，瞬间降到 80 MB，这体验让我彻底信了稀疏矩阵。

为什么空间单细胞数据天然适合稀疏表示

想象一下你手里那张空间表达矩阵：行是基因，列是 spot。绝大多数基因只在少数 spot 里有表达，尤其是那些细胞类型标记基因。你想想看，一个 spot 里能同时高表达的基因通常不超过 1000 个，而总基因数可能是 20000。这意味着矩阵的稀疏度（零元素占比）普遍在 95% 以上。

• 直接存稠密矩阵：20000 × 5000 × 4 bytes = 400 MB（float32）
• 存为 CSR/CSC：非零元素约 5000×1000=5M 个，加上行偏移和列索引，大约 5M×(4+4) + 5001×4 ≈ 40 MB
• 节省 90% 内存，而且索引结构本身还能加速矩阵运算

所以说，处理空间单细胞数据的第一个习惯，就是别用 Numpy 的 ndarray 去存原始计数。我早期犯过这个错，结果 swap 分区疯狂读写，差点把工作站搞死。

CSR 与 CSC：两种最常用的稀疏格式

scipy.sparse 里 CSR（Compressed Sparse Row）和 CSC（Compressed Sparse Column）是对称的兄弟。你只需要搞明白一个，另一个自然就通了。

对比维度	CSR	CSC
存储结构	按行压缩，存行偏移 + 列索引 + 数据	按列压缩，存列偏移 + 行索引 + 数据
切行效率	O(非零元素数)，极快	慢，需要逐列搜索
切列效率	慢	O(非零元素数)，极快
矩阵乘法（左乘）	高效	低效
矩阵乘法（右乘）	低效	高效
典型使用场景	基因×spot（行是基因）	spot×基因（行是 spot）

我个人习惯：如果后续要频繁按基因取行（比如筛选高变基因），就用 CSR；如果频繁按 spot 取列（比如提取某个 spot 的表达谱），就用 CSC。在 Scanpy 的底层实现里，anndata 的 X 矩阵默认是 CSR——因为大部分操作是按基因（obs）来过滤的。

实战：用 scipy.sparse 构建和操作

我不会一上来就讲 scipy.sparse 的全部 API，那太枯燥了。咱们直接拿一个模拟的空间表达矩阵来走一遍流程。

import numpy as np
from scipy import sparse

# 模拟一个 20000×5000 的稀疏矩阵，密度 5%
np.random.seed(42)
n_genes = 20000
n_spots = 5000
density = 0.05
n_nonzero = int(n_genes * n_spots * density)

# 随机生成非零坐标和值
row = np.random.randint(0, n_genes, n_nonzero)
col = np.random.randint(0, n_spots, n_nonzero)
data = np.random.poisson(lam=5, size=n_nonzero).astype(np.float32)

# 构造 CSR 矩阵
csr_mat = sparse.csr_matrix((data, (row, col)), shape=(n_genes, n_spots))
print(f"CSR 内存占用: {csr_mat.data.nbytes + csr_mat.indices.nbytes + csr_mat.indptr.nbytes} bytes")
print(f"非零元素数量: {csr_mat.nnz}")

# 构造 CSC 矩阵（从 CSR 转换）
csc_mat = csr_mat.tocsc()
print(f"CSC 内存占用: {csc_mat.data.nbytes + csc_mat.indices.nbytes + csc_mat.indptr.nbytes} bytes")

输出结果会让你直观看到内存差距。注意这里我用了 float32，而不是默认的 float64——这也是个优化点，稍后会提。

接下来演示常见的行/列操作：

# 取第 0 个基因的表达（CSR 直接切行）
gene_0 = csr_mat[0, :]  # 返回的还是一个 CSR 矩阵
# 取第 100 个 spot 的表达（CSR 切列要转成 CSC）
spot_100 = csc_mat[:, 100]  # 返回 CSC，快

# 矩阵乘法：计算基因共表达（左乘）
gene_corr = csr_mat @ csr_mat.T  # 结果形状 20000×20000，小心内存
# 实际不要这样干，用稀疏的 sparse.csgraph 或者只算 top-k

# 行求和：每条基因的表达总量
gene_sum = np.array(csr_mat.sum(axis=1)).flatten()
print(f"前 10 个基因表达量: {gene_sum[:10]}")

内存对比实测（我踩过的坑）

我用自己项目中的真实数据跑过一次对比：10× Visium 小鼠脑切片，基因数 19,480，spot 数 2,836。

表示格式	内存占用	稀疏度	备注
稠密 float64	19,480×2,836×8 ≈ 442 MB	0%	完全不可接受
稠密 float32	221 MB	0%	还是太大
CSR float64	约 38 MB	92.3%	可用
CSR float32	约 22 MB	92.3%	推荐方案

你看，从 442 MB 降到 22 MB，节省了 20 倍。而且最关键的是，CSR 格式下做 sum、mean、slice 都不会生成中间稠密矩阵。我当初用稠密矩阵跑 Louvain 聚类，直接卡了 3 分钟；换 CSR 后 15 秒搞定。嗯，这只是格式带来的收益，还没算并行。

SVG 图：稀疏矩阵表示知识体系

一些实用建议与避坑

说了这么多，最后给出几条我在项目里总结的黄金法则：

• 数据存为 CSR，计算时按需转格式。比如你需要做按列操作，就用 mat.tocsc() 转一下，但注意转换本身有开销，别在循环里反复转。
• 永远用 float32。空间单细胞的表达值基本都是整数计数，float32 的精度完全够用。我曾经为了省事用默认 float64，结果内存翻倍，速度也慢了一半。
• 避免稀疏矩阵直接参与广播。例如 csr_mat + np.array([1,2,3]) 会触发转换为稠密。应该手动对齐维度，或者用 sparse 的 broadcast_to 工具。
• 用 .nnz 属性监控稀疏度。如果稀疏度降到 50% 以下，可能稠密格式反而是更好的选择（但空间单细胞几乎不会出现这种情况）。