数据准备：空间转录组数据的读取与预处理

把原始数据喂进 Visium 聚类 pipeline 之前，得先过一遍「数据准备」这关。

说白了，这一步干的就三件事：读进来、洗干净、存好档。别小看这步，我见过太多项目栽在数据预处理上——基因表达矩阵里混了一堆死细胞、空液滴的噪音，后面聚类结果怎么看怎么不对。

嗯，咱们今天就把这块啃透。

读取 Visium 标准输出

10x Genomics 的 Visium 平台会把测序结果打包成几个固定文件。

• filtered_feature_bc_matrix.h5 —— 过滤后的表达矩阵
• raw_feature_bc_matrix.h5 —— 原始表达矩阵
• spatial/tissue_lowres_image.png —— 组织切片低分辨率图
• spatial/tissue_positions_list.csv —— 斑点坐标表
• spatial/scalefactors_json.json —— 缩放因子等元信息

我个人习惯用 Scanpy 做主力加载工具，它内置了 Visium 数据读取接口，不用手写 I/O 循环。

import scanpy as sc
import anndata

# 读取 Visium 输出目录
adata = sc.read_visium(
    path='./data/visium_sample/',
    count_file='filtered_feature_bc_matrix.h5',
    load_images=True
)

print(adata)
# AnnData object with n_obs × n_vars = 4992 × 33538
#    obs: 'in_tissue', 'array_row', 'array_col'
#    var: 'gene_ids', 'feature_types', 'genome'
#    uns: 'spatial', 'spatial_neighbors'
#    obsm: 'spatial'

读完后，adata 对象里会多出 obsm['spatial'] 这一项，存的就是每个斑点的 X/Y 像素坐标。我建议你立即检查一下这张表，确认坐标数量和窗口观察数一致。

质量过滤的 SOP

空间转录组数据里埋了不少坑。最常见的三个雷：

1. 低质量斑点 —— 捕获的基因太少，基本是空液滴或死细胞
2. 线粒体基因占比过高 —— 细胞膜破损，线粒体内容物泄漏
3. 双细胞污染 —— 一个斑点里捕获了两个以上细胞

我是按下面这个流程筛的，你直接套用就行。

# 1. 计算 QC 指标
adata.var['mt'] = adata.var_names.str.startswith('MT-')
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False)

# 2. 查看分布
sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],
             jitter=0.4, multi_panel=True)

# 3. 应用过滤条件
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)

# 手动过滤双细胞和线粒体
adata = adata[adata.obs['n_genes_by_counts'] < 6000, :]
adata = adata[adata.obs['pct_counts_mt'] < 20, :]

print(f'剩余斑点数: {adata.n_obs}')
# 剩余斑点数: 4587

你想想看，如果不做这步过滤，下游聚类会把噪声当成信号。我踩过这个坑——有一回偷懒没筛 mt%，结果聚类算法分出了一整群「线粒体富集」的假簇，差点带着整个解读方向跑偏。

基因与 UMI 的二次筛选

过滤完斑点，还得回头清理基因维度。

有些基因在整个组织里只有零星几个 UMI 计数，它们对聚类贡献极低，却会增加计算负担。我的经验是：至少要在 3 个斑点里有非零表达才保留。

# 保留至少在 3 个细胞中表达的基因
sc.pp.filter_genes(adata, min_cells=3)

# 进一步过滤：表达量极低的基因也删掉
# 计算每个基因在所有斑点的平均表达量
gene_means = adata.X.mean(axis=0).A1  # 稀疏矩阵转 array
genes_to_keep = gene_means > 0.01
adata = adata[:, genes_to_keep]

print(f'剩余基因数: {adata.n_vars}')
# 剩余基因数: 14220

归一化与对数变换

原始 UMI 计数受测序深度、斑点大小的影响很大。

不归一化直接做 PCA，第一主成分大概率就是测序深度，那聚类就完蛋了。所以一定要做两步标准操作。

# 1. 全局尺度标准化：每个斑点的总计数归一化到中位数
sc.pp.normalize_total(adata, target_sum=1e4)

# 2. 对数变换（ln(1 + x)）
sc.pp.log1p(adata)

你可能会问，为什么非要用 ln(1+x)？因为原始计数是偏移很大的泊松分布，对数变换能让高表达基因和低表达基因的方差更均匀，PCA 和聚类算法才能正常工作。

高变基因（HVG）选择

空间转录组的基因数通常在上万级别，全量跑 PCA 很慢，而且噪声太多。

我的做法是挑选 2000~5000 个「高变基因」——说白了就是在不同斑点间表达量波动最大的那些基因。这些基因才是区分细胞类型的主力。

# 挑选前 3000 个高变基因
sc.pp.highly_variable_genes(adata, n_top_genes=3000, flavor='seurat_v3')

# 如果 flavor='seurat_v3' 报错，试一下默认的 'seurat'
sc.pp.highly_variable_genes(adata, n_top_genes=3000, flavor='seurat')

# 保存高变基因标记，后续 PCA 只跑这些基因
adata.var['highly_variable'] = adata.var['highly_variable']
adata_hvg = adata[:, adata.var['highly_variable']].copy()

数据准备管线总览

把上面几步串起来，整个数据准备流程就是下面这张图。你可以把它贴在项目文档里当 SOP。

按这个流程走下来，你的 adata 对象就准备好接收 PCA、聚类参数调优了。