📚 课程目录第4章空间基因表达矩阵构建

空间基因表达矩阵构建：从稀疏矩阵到工具实战

好，这一章我们聊聊基因表达矩阵。说白了，就是把每个空间点的基因计数整理成一个规整的数据结构。听起来简单？我刚开始做空间转录组的时候，在这步栽过不少跟头。你想想看，一个10x Visium芯片就有5000个spot，每个spot检测2万个基因，算下来就是1亿个数值。这还没算更高分辨率的平台。

嗯，这里要注意——绝大多数的数值都是零。为什么？因为单个spot里只捕获了几十到几百个细胞，大部分基因压根没表达。这就是我们要用稀疏矩阵的核心原因。

稀疏矩阵：空间转录组的基石

常规的稠密矩阵，存1亿个浮点数得占多少内存？按float32算，大约400MB。单个样本还能忍，但空间转录组动辄几十个样本，再叠加下游降维、聚类、空间邻域计算……内存很快就撑爆了。

稀疏矩阵只存储非零元素，以及它们的位置信息。常见的有三种格式：

COO（坐标列表）：存(row, col, value)三元组。直观，适合增量构建。
CSR（压缩稀疏行）：行指针+列索引+数值。矩阵乘法快，是scipy默认格式。
CSC（压缩稀疏列）：列指针+行索引+数值。列切片效率高。

个人经验：我习惯用CSR格式做下游分析。大部分空间转录组的操作是按spot（行）来处理的，比如找邻居、归一化。CSR对行操作极度友好。

来看看实际代码。假设你已经有了一个counts字典：

import scipy.sparse as sp
import numpy as np

# 模拟数据：5000个spot，20000个基因
n_spot, n_gene = 5000, 20000
# 假设只有5%的位置有表达
density = 0.05
n_nonzero = int(n_spot * n_gene * density)

# 随机生成非零位置和计数
row = np.random.randint(0, n_spot, n_nonzero)
col = np.random.randint(0, n_gene, n_nonzero)
data = np.random.poisson(3, n_nonzero)

# 构建COO矩阵
coo = sp.coo_matrix((data, (row, col)), shape=(n_spot, n_gene))

# 转成CSR，这是最常用的格式
csr = coo.tocsr()

# 看看内存节省了多少？
import sys
dense_memory = sys.getsizeof(np.zeros((n_spot, n_gene), dtype=np.float32))
sparse_memory = sys.getsizeof(csr.data) + sys.getsizeof(csr.indices) + sys.getsizeof(csr.indptr)

print(f"稠密矩阵: {dense_memory / 1024 / 1024:.2f} MB")
print(f"稀疏矩阵: {sparse_memory / 1024 / 1024:.2f} MB")

小提示：上面的代码只是演示原理。真实项目里，你一般不会手动构建稀疏矩阵。主流工具已经封装好了。

Seurat vs Squidpy：工具选型

这个选择题我当年也纠结过。后来两个都用熟了，发现没有绝对的优劣，关键看你的数据规模和分析偏好。

维度	Seurat (R语言)	Squidpy (Python)
数据对象	SeuratObject + SpatialImage	AnnData
稀疏矩阵存储	dgCMatrix (内置)	scipy.sparse.csr_matrix
空间邻接计算	基于距离 + 组织学图像	基于距离 + 组织学 + 图论
与深度学习对接	需通过 h5Seurat 或 rpy2 桥接	原生对接 PyTorch/TensorFlow
社区成熟度	极高，教程丰富	快速增长，论文标配

我个人的建议是：

如果团队以R为主，分析流程偏经典——选Seurat。它把空间数据整合得相当顺手。
如果要对接深度学习模型，或者做自定义分析——选Squidpy。Python生态实在太香了。

Seurat实战：从空间数据到SeuratObject

Seurat从4.0版本开始内置了空间数据支持。来看看标准流程：

# R语言示例
library(Seurat)
library(SeuratObject)

# 读取空间数据
# 假设你的数据是10x Visium格式
visium_data <- Read10X_Spatial(
  data.dir = "filtered_feature_bc_matrix/",  # 表达矩阵
  images = "spatial/",                         # 组织学图像
  assay = "Spatial"                            # 指定assay名
)

# 创建Seurat对象
seurat_obj <- CreateSeuratObject(
  counts = visium_data$Spatial,
  assay = "Spatial"
)

# 添加空间坐标
seurat_obj <- AddTissueCoordinates(
  object = seurat_obj,
  image = visium_data$image,
  scale.factor = 60  # 10x Visium默认缩放因子
)

# 查看内部存储：本质上是dgCMatrix
class(GetAssayData(seurat_obj, slot = "counts"))
# [1] "dgCMatrix"

曾经踩过的坑：AddTissueCoordinates这一步很容易出错。空间坐标的单位是像素，不是微米。如果你后续要做距离计算，记得先转换。我有个项目跑了三天才发现距离全算错了，就是因为没处理这个单位问题。

Squidpy实战：AnnData + 稀疏矩阵

Squidpy基于AnnData，而AnnData的核心矩阵就是CSR。来看看怎么构建：

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

# 读取空间数据（假设是Visium H5格式）
adata = sc.read_visium("spatial/", library_id="sample1")

# 检查表达矩阵的存储格式
print(type(adata.X))
# 

# 看看稀疏比例
n_total = adata.shape[0] * adata.shape[1]
n_nonzero = adata.X.nnz
print(f"稀疏度: {n_nonzero / n_total * 100:.2f}%")

# 添加空间邻域关系
sq.gr.spatial_neighbors(adata, coord_type="grid", n_neighs=6)

# 查看邻接矩阵（也是稀疏的！）
print(type(adata.obsp["spatial_distances"]))
# 

# 计算空间变量基因
sq.gr.spatial_autocorr(adata, mode="moran", n_perm=1000)

我个人习惯用AnnData格式。为什么？因为它在设计时就把“矩阵+元数据”当成一个整体，不会像Seurat那样经常需要手动对齐索引。而且，AnnData的.obsm和.obsp可以直接存PCA、邻接矩阵这些中间结果，非常方便。

矩阵构建的避坑指南

我在项目中遇到过几个经典问题，写在这里供你参考：

基因名重复：不同批次的数据合并时，同一基因可能有别名（比如"ENSG000001..."和"TP53"）。先做基因注释统一，再建矩阵。
Spot过滤不及时：低质量的spot（总counts < 100）会拖慢分析，还会引入噪音。建议在构建矩阵时就过滤掉。
内存泄漏：如果你在循环里反复构建新的稀疏矩阵，记得手动gc.collect()。Python的垃圾回收有时候不太及时。

一个小技巧：构建完矩阵后，立刻做一次完整性检查——确保矩阵的稀疏模式、行名、列名都与预期一致。我后来把检查写成了pytest用例，每次跑数据前自动验证。

数据质量控制与矩阵修剪

矩阵构建完，不等于可以直接丢进模型。质量过滤是必要的一步：

# 基于统计的过滤策略
import numpy as np
import scipy.sparse as sp

def filter_matrix(csr_matrix, min_genes=200, max_genes=5000, min_counts=500):
    """
    过滤低质量spot和低表达基因
    
    参数:
        csr_matrix: spot × gene 的CSR矩阵
        min_genes, max_genes: 每个spot检测到的基因数范围
        min_counts: 每个spot的最小UMI计数
    """
    # 按spot过滤
    n_genes_per_spot = np.diff(csr_matrix.indptr)
    n_counts_per_spot = np.array(csr_matrix.sum(axis=1)).flatten()
    
    valid_spots = (
        (n_genes_per_spot >= min_genes) & 
        (n_genes_per_spot <= max_genes) & 
        (n_counts_per_spot >= min_counts)
    )
    csr_filtered = csr_matrix[valid_spots, :]
    
    # 按基因过滤（至少在3个spot中有表达）
    gene_counts = np.diff(csr_filtered.tocsc().indptr)
    valid_genes = gene_counts >= 3
    csr_filtered = csr_filtered[:, valid_genes]
    
    return csr_filtered, valid_spots, valid_genes

嗯，这里要注意过滤的顺序。我建议先过滤spot，再过滤基因。因为无效的spot会给基因统计带来噪音，反过来会让好基因也被误杀。

最后，来看看这一章的知识结构。我画了一张流程图，帮你理清从原始数据到标准分析对象的完整链路：

从上图可以看出，核心路径只有几步：原始数据 → 稀疏矩阵 → 工具适配 → 质量控制 → 标准化 → 可分析对象。每一步都有对应的稀疏矩阵操作，理解了这些，你就掌握了空间转录组数据处理的基础。