企业官网建设流程全解析

标题写成 "10 x visum ST 预处理"，玩过空间转录组的朋友一眼就能看出来，这是 10x Visium ST 的笔误或者说简写，ST 就是 Spatial Transcriptomics。我最近刚好完整处理了一批 Visium ST 的原始数据，从读入、预处理到 Python 可视化收尾，把每一步的参数选择、代码逻辑和踩过的坑都整理了一遍。这篇文章就是这份实操记录，目标读者是两类人：刚接触空间转录组、准备用 Python 跑 Visium ST 预处理和可视化的新手；以及已经能跑通 Scanpy 流程、但想搞清楚每一步为什么这么做、遇到问题怎么排查的进阶用户。

1. 项目拆解：Visium ST 预处理到底要做什么

1.1 一份 Visium ST 数据里有什么

拿到 Visium ST 数据之后，先别急着写代码。找一个样本的输出目录，把文件结构从头到尾看一遍，这是最容易被跳过、但恰恰最值得花时间的步骤。一个标准的 Visium ST 输出目录大概长这样：

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

其中最关键的有两块：表达矩阵和空间信息。filtered_feature_bc_matrix.h5 是过滤后的表达矩阵，barcode 对应组织切片上的每一个 spot（捕获点），feature 就是基因。而 spatial 目录下那一堆文件，才是 Visium ST 跟普通单细胞数据最不一样的地方。

tissue_positions_list.csv 里记录着每个 spot 的条形码和坐标。新手最容易在这里犯晕，因为里面有两套坐标系。array_row 和 array_col 是捕获芯片上的阵列坐标，走的是六边形网格布局，跟组织图像没有直接关系；pxl_col_in_fullres 和 pxl_row_in_fullres 才是 spot 在全分辨率组织图像上的像素坐标。做可视化对齐时，真正用到的是后者。

scalefactors_json.json 记录了一组缩放因子。全分辨率图像太大，不适合直接展示和计算，所以官方流程提供了高分辨率和低分辨率两张组织图。要把全分辨率像素坐标换算到低分辨率图上，就必须用缩放因子做一次坐标变换。读入数据时如果忽略这一步，画出来的 spot 位置一定是偏的，甚至整片飞出组织轮廓。

1.2 预处理流程的全貌与定位

Visium ST 的预处理，本质上是一条跟单细胞分析高度相似的流水线：质控、过滤、标准化、高变基因筛选、PCA、邻域图构建、聚类，最后是可视化。但有两个关键差异，决定了你不能无脑照搬单细胞教程。

第一，质控的单位不同。单细胞质控的对象是细胞，Visium ST 质控的对象是 spot。一个 spot 可能捕获了多个细胞的 RNA，也可能是组织边缘或空白区域，几乎没捕获到东西。这意味着 spot 的基因数和 UMI 数分布会比单细胞更宽，不能直接套用"500 到 5000 个基因"那套经验阈值。

第二，多了一条图像对齐的链路。分析结果最终要映射回组织切片上，涉及坐标变换和图像处理，这在单细胞分析里完全不存在。再加上空间转录组技术本身还在快速迭代，不同批次、不同样本的切片质量和 spot 密度差异很大，预处理阶段就要把样本质量差异识别出来，避免下游分析被个别坏样本带偏。

所以我的建议是：预处理阶段不要把目标定得太高，先把四件事跑通——数据正确读入、质控合理、聚类结果可解释、空间可视化对得上。下游的细胞通讯、空间共定位、区域差异分析都是后话，前提是前面的数据质量过关。下面进入实操。

2. 环境准备与数据读入

2.1 工具链选型与安装

做 Visium ST 预处理的 Python 工具不少，但主力还是 Scanpy。原因很简单：它内置了 read_visium 函数，把 h5 表达矩阵、空间坐标、图像信息一次性读入 AnnData 对象，省掉大量手写解析代码的功夫。配合 Squidpy 做空间统计和邻域分析，Matplotlib 做图像输出，基本能覆盖全流程。

环境建议用 conda 或 mamba 创建：

conda create -n visium python=3.9 -y conda activate visium pip install scanpy squidpy matplotlib seaborn h5py

Scanpy 版本建议用 1.9 或更新版本，旧版本对 read_visium 的兼容性和新版 h5 文件格式支持都不太好。这里补一句：如果机器内存不大，把 anndata 和 scanpy 都升级到最新版，旧版 anndata 在读写大文件时会有明显的内存开销，这是我在处理 8G 内存笔记本时踩过的实坑。

2.2 read_visium 读入与手动解析

最省事的就是 Scanpy 自带接口：

import scanpy as sc import pandas as pd import numpy as np import matplotlib.pyplot as plt import json adata = sc.read_visium( path="path/to/outs/", count_file="filtered_feature_bc_matrix.h5", library_id="sample1", source_image_path=None, )

这里有几个参数值得讲清楚。path 指向的是包含 spatial 目录的文件夹；library_id 是自己给样本起的名字，多个样本时可以用不同 ID 区分；source_image_path 通常保持 None，让 Scanpy 自己去 spatial 目录里找图像。读入之后立刻检查结构：

print(adata) print(adata.obs.head()) print(adata.uns["spatial"])

一切正常的话，adata.obs 里会看到 array_row、array_col、pxl_col_in_fullres、pxl_row_in_fullres 等列，adata.uns["spatial"] 里保存着图像和缩放因子信息。

提示：如果读入时报错或者 obs 里没有坐标列，先确认输出目录的版本。较老的 10x 格式里坐标文件名是 tissue_positions_list.csv 且没有表头，新格式可能是 tissue_positions.csv 且带表头。Scanpy 会自动识别，但你手动解析时必须分清，否则列名会整体错位。

另一种方式是直接用 h5py 手动读 h5 文件，虽然麻烦，但能帮你理解数据本质：

import h5py with h5py.File("path/to/outs/filtered_feature_bc_matrix.h5", "r") as f: print(list(f.keys())) grp = f["matrix"] barcodes = grp["barcodes"][:] data = grp["data"][:] indices = grp["indices"][:] indptr = grp["indptr"][:] shape = grp["shape"][:]

这就是稀疏矩阵的 CSC 格式。barcodes 对应 spot，features 在另一个 group 里存着基因名和类型。我建议新手至少手动读一次，理解了这个结构，将来遇到 Scanpy 版本不兼容，或者想读取 raw_feature_bc_matrix.h5 做自定义过滤时，就不会束手无策。坐标文件同理：

positions = pd.read_csv( "path/to/outs/spatial/tissue_positions_list.csv", header=None, names=["barcode", "in_tissue", "array_row", "array_col", "pxl_col_in_fullres", "pxl_row_in_fullres"], )

注意这里的 header=None 是针对旧版无表头文件。如果文件自带表头，就要用 header=0，并按实际列名映射。两套文件读完后，可以用 barcode 把表达矩阵和坐标 merge 起来，这就是 read_visium 在背后做的事。

3. 质控过滤与标准化实战

3.1 质控指标：从单细胞思维切换到 spot 思维

质控的第一步是算 QC 指标，然后判断哪些 spot 质量差。核心指标有三个：总 UMI 数（total_counts）、检测到的基因数（n_genes_by_counts）、线粒体基因比例（pct_counts_mt）。

adata.var["mt"] = adata.var_names.str.startswith("MT-") sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], percent_top=None, inplace=True)

很多从单细胞转过来的人会疑惑：组织切片上的 spot 又不是完整细胞，线粒体比例还有意义吗？答案是：有，但解读方式不同。如果某个 spot 的线粒体比例异常高，通常意味着这个 spot 捕获到的 RNA 质量差、细胞破碎严重，或者正处于坏死区域。它不像单细胞那样用来筛掉"垂死细胞"，而是用来标记异常区域，辅助判断组织质量。

画图看分布是关键一步：

sc.pl.violin(adata, ["total_counts", "n_genes_by_counts", "pct_counts_mt"], multi_panel=True, jitter=0.3)

也可以直接看散点分布：用 total_counts 和 n_genes_by_counts 做双变量散点，能直观看到是否有明显的离群 spot 群。对照 tissue_positions_list.csv 里的 in_tissue 字段，你会发现大量低计数的 spot 基本都落在组织外或边缘。这也是 Visium 数据质控的关键技巧：in_tissue=0 的 spot 默认就不该参与分析。

3.2 过滤阈值怎么定

过滤没有绝对标准，但有三条经验路径可以组合使用。

路径一：绝对阈值。比如认为一个 spot 至少要有 200 个基因、500 个 UMI，基因至少要在 5 个 spot 中有表达，才值得保留：

sc.pp.filter_cells(adata, min_genes=200) sc.pp.filter_cells(adata, min_counts=500) sc.pp.filter_genes(adata, min_cells=5)

路径二：百分位数。如果样本整体质量不错，取分布的低百分位作为下限更稳妥。用 1% 或 5% 分位数做下限，能温和地去掉尾部低质量 spot，同时不误伤组织边缘的真实信号：

lower_count = np.percentile(adata.obs["total_counts"], 1) lower_gene = np.percentile(adata.obs["n_genes_by_counts"], 1) adata = adata[adata.obs["total_counts"] > lower_count, :] adata = adata[adata.obs["n_genes_by_counts"] > lower_gene, :]

路径三：MAD 法。用中位数加减 3 倍绝对中位差来判定离群点。这种方法对"整体分布偏移"更稳健，适合在不同批次间做标准化处理，但用在单样本 Visium 数据上时，注意不要在高质量样本里把阈值卡得太严。

注意：过滤阈值定了之后，一定要用空间可视化去验证。把 total_counts 画到组织切片上，如果低计数 spot 集中在组织边缘或空白区域，说明过滤方向是对的；如果组织内部出现大片低计数区域，那可能是切片本身有气泡或组织坏死，此时不要盲目提高阈值，而应该先跟实验记录对照。

3.3 标准化、对数化与高变基因

过滤完成后进入标准化环节：

sc.pp.normalize_total(adata, target_sum=1e4) sc.pp.log1p(adata) sc.pp.highly_variable_genes(adata, n_top_genes=2000, flavor="seurat")

normalize_total 把每个 spot 的计数缩放到总和 10000，相当于"每 spot 测了 1 万 UMI 时各基因会是多少"，消除测序深度差异。target_sum 选 1e4 是惯例，实际操作中也可以根据总 UMI 中位数调整，但不建议随便改，因为后续差异表达和可视化的很多阈值习惯都建立在 1e4 这个尺度上。

log1p 是 ln(x+1) 变换，压缩动态范围，让强表达基因不至于在后续 PCA 中完全压过弱表达基因。这里需要记住顺序：如果要用 seurat_v3 这个 flavor 选高变基因，它要求输入原始计数，而非对数化数据。我一般这样处理：

adata.layers["counts"] = adata.X.copy() sc.pp.normalize_total(adata, target_sum=1e4) sc.pp.log1p(adata) sc.pp.highly_variable_genes(adata, layer="counts", flavor="seurat_v3", n_top_genes=2000)

把原始计数备份到 layer 里，再用 layer 参数调用 seurat_v3，两全其美。选完高变基因后，通常会把数据裁剪到高变基因子集，再做标准化：

adata = adata[:, adata.var["highly_variable"]].copy() sc.pp.scale(adata, max_value=10)

scale 就是 z-score 标准化，让每个基因均值为 0、方差为 1，避免高表达基因在降维时占据主导。max_value=10 把极端值截断，这个细节很多人会忽略，但它在后续 PCA 时能显著减少离群 spot 对主成分方向的影响。

4. 降维与聚类：空间数据的邻域效应

4.1 PCA 与邻域图构建

高变基因选完、数据缩放之后，进入降维和聚类环节：

sc.tl.pca(adata, svd_solver="arpack", n_comps=30) sc.pp.neighbors(adata, n_neighbors=10, n_pcs=30) sc.tl.umap(adata, min_dist=0.3, random_state=0) sc.tl.leiden(adata, resolution=0.5, random_state=0)

PCA 这一步用 arpack 求解器就够了，Visium 数据通常几千个 spot、2000 个高变基因，样本量不大，不需要 full 或 randomized 求解器。n_comps 选 30 是折中方案：太少会丢失信息，太多会引入噪声。跑完可以看一眼方差贡献曲线：

sc.pl.pca_variance_ratio(adata, n_pcs=50)

通常前 20 到 30 个主成分就能解释大部分方差，如果曲线在第 10 个主成分附近就进入平台期，可以适当把 n_pcs 降到 15 左右。

neighbors 是聚类和 UMAP 共用的"图构建"步骤。它基于 PCA 空间计算 spot 间的相似度，构建 k 近邻图。n_neighbors 默认 10 到 15 都行，数值越小越能捕捉局部结构、但噪声也越大；数值越大图越平滑、但会模糊掉细小的空间区域。空间转录组数据的特殊性在于，spot 之间天然存在空间邻域关系，所以很多人会在 neighbors 之后额外叠加空间距离约束，比如用 Squidpy 构建空间邻域图，或者直接用空间聚类方法。但作为基线流程，Leiden 聚类已经足够先看整体结构。

4.2 分辨率参数与聚类结果解读

Leiden 的 resolution 参数直接控制聚类的粗细。0.5 是通用默认值，但对空间转录组数据，我通常会从 0.3 开始尝试，逐步提高到 1.0，观察聚类结果是否与组织解剖结构对应。比如脑切片，理想的聚类结果应该能大致区分皮层、白质、海马等区域。

有个实用技巧：用一个 for 循环依次生成不同 resolution 的 Leiden 聚类，然后批量画空间散点图对比，比反复手动调整效率高得多。只看数字和 marker 基因很难判断"这个聚类好不好"，但在空间图上，聚类边界是否跟组织解剖结构吻合、有没有出现大量"斑点状"的碎片聚类，一眼就能看出来。

UMAP 在这条流程里的角色要摆正。它用于展示全局转录组相似性，帮助发现"哪几个聚类在表达谱上更接近"，但它不能替代空间图的位置信息。我在实际报告中，一般把 UMAP 和空间散点图并排展示：UMAP 反映转录组状态，空间图反映组织位置，两者结合才能讲完整的生物学故事。

提示：每次跑 UMAP 和 Leiden 都要固定 random_state，否则每次运行结果都会有细微差别，这在多人协作或者复现分析时很要命。我见过不止一次，因为没固定随机种子，第二天重新跑一遍流程，聚类编号全变了，下游所有结果都要返工。

5. Python 可视化实战

5.1 空间坐标与图像对齐的核心逻辑

很多人卡在"为什么我的 spot 和图像对不上"这个问题上。理解了对齐逻辑，问题就解决了一半。

Visium 的原始坐标是 spot 在全分辨率图像上的像素位置（pxl_col_in_fullres、pxl_row_in_fullres）。全分辨率图像动辄几万像素，不适合直接展示，所以官方提供了降采样后的 hires 和 lowres 图像，并给出缩放因子。要把 spot 画到 hires 图上，实际做的是这样的变换：

with open("path/to/outs/spatial/scalefactors_json.json") as f: scalefactors = json.load(f) positions = pd.read_csv("path/to/outs/spatial/tissue_positions_list.csv", header=None) scalef = scalefactors["tissue_hires_scalef"] x = positions[4].values * scalef y = positions[5].values * scalef img = plt.imread("path/to/outs/spatial/tissue_hires_image.png") fig, ax = plt.subplots(figsize=(8, 8)) ax.imshow(img) ax.scatter(x, y, s=1, c="red", alpha=0.3)

注意 x 用的是 pxl_col（列坐标），y 用的是 pxl_row（行坐标），不要搞反。图像坐标系里，row 是 y、col 是 x，直接用就行。spot 的直径也要缩放，在图上画圆点大小时可以参考 spot_diameter_fullres 乘以缩放因子，一般是 50 到 100 像素左右。如果只关注组织内的 spot，可以把 in_tissue 列筛出来再画，图面会干净很多。

Scanpy 的 read_visium 已经在内部帮你做了这套变换，并且把 coordinates 存进了 adata.obsm["spatial"]。动手实现一遍的唯一目的是：当结果对不上时，你能自己排查是哪一步出了问题。

5.2 核心空间可视化命令

有了 AnnData 对象，最常用的空间可视化是一行命令：

sc.pl.spatial( adata, color="leiden", img_key="hires", alpha=0.6, size=1.2, frameon=False, title="Spatial clusters", )

这里几个参数值得细说。img_key 选 hires 还是 lowres，取决于展示需求：hires 组织细节丰富，适合放到论文图里；lowres 文件小、加载快，适合快速探索。alpha 是点透明度，控制在 0.5 到 0.7 之间最合适，太透明看不清聚类边界，太不透明会挡住组织形态。size 是 spot 的大小，1.0 到 1.5 是常用区间，spot 太大会互相重叠，太小又看不出区域轮廓。

质控指标同样可以画到空间图上，这是验证预处理效果的关键动作：

sc.pl.spatial(adata, color="total_counts", img_key="hires", size=1.2) sc.pl.spatial(adata, color="n_genes_by_counts", img_key="hires", size=1.2) sc.pl.spatial(adata, color="pct_counts_mt", img_key="hires", size=1.2)

如果 total_counts 空间分布跟组织形态高度吻合，低计数区域集中在背景或边缘，说明样本质量正常、过滤处理得当。如果组织内部出现明显的低计数空洞，就要警惕切片质量问题。

自定义空间图也很简单。拿到 obsm["spatial"] 里的坐标，直接用 Matplotlib 的 scatter 就能画：

coords = adata.obsm["spatial"] plt.figure(figsize=(8, 8)) plt.scatter(coords[:, 0], coords[:, 1], c=adata.obs["leiden"].astype(str), s=8, alpha=0.7)

这种自定义方式灵活，比如只画某个聚类、按基因表达值着色、或者叠加不同样本的边界，都是 sc.pl.spatial 不太方便做的。

5.3 Marker 基因的空间分布

聚类完成后，通常需要找每个聚类的 marker 基因，然后看它们的空间分布是否解释得通：

sc.tl.rank_genes_groups(adata, groupby="leiden", method="wilcoxon") sc.pl.rank_genes_groups(adata, n_genes=8, sharey=False)

wilcoxon 秩和检验是常用的 marker 检测方法，它对表达分布的非正态性有较好的鲁棒性。跑完之后，把感兴趣的基因直接画到空间图上：

sc.pl.spatial( adata, color=["Mbp", "C1qa", "Sst"], img_key="hires", ncols=3, size=1.2, frameon=False, )

sc.pl.spatial 的 color 参数可以直接传基因名，它会把基因表达量映射到颜色梯度。看 marker 基因的空间分布，是验证聚类生物学意义最直观的一步。比如脑组织里，髓鞘基因富集的聚类应该落在白质区域，小胶质细胞 marker 富集的聚类应该分布在整个脑区，如果某个已知的区域 marker 出现在完全错误的位置，那聚类参数或者数据质量可能有问题。

注意：画基因表达时，如果基因在数据里过滤掉了，会直接报错。用 adata.var_names 确认基因名存在，注意基因名大小写和版本（比如线粒体基因是 MT- 前缀还是 mt- 前缀），不同参考基因组版本会影响这个。

6. 常见问题与排查技巧实录

6.1 典型问题速查表

我在处理 Visium ST 数据时，把最容易遇到的问题和排查思路整理成了下面的表格，基本都是真实踩过的坑：

6.2 独家避坑技巧

除了表格里的排查思路，还有几个小技巧值得单独拿出来说。

第一，每完成一步预处理就保存一次 AnnData。我习惯在过滤后、标准化后、聚类后分别保存 h5ad 文件，文件命名带步骤标记。Visium 数据动辄几个 GB，重新跑一遍过滤要花不少时间。更重要的是，下游分析经常要回到中间状态重新试参数，有保存点就不用从头再来。

第二，先用一个小样本或者裁剪区域跑通全流程。如果同时处理多个样本，不要一上来就跑全部数据。抽一个代表性样本，把 read_visium 到最终可视化的全流程走通，确认参数合理后再批量应用。批量处理时写一个 for 循环，把 library_id 和输出路径参数化，可以在半小时内完成几十个样本的预处理。

第三，注意多个切片样本的批次效应。Visium 实验常涉及多个切片、多个个体，即使来自同一组织，不同切片间的测序深度和 spot 捕获效率也可能有差异。预处理阶段每个样本单独跑，得到聚类后如果想做跨样本比较，需要做批次整合或者至少用多个样本的联合高变基因重新分析。不能把不同样本的原始表达矩阵直接拼接后跑聚类，那会把批次差异和生物学差异混在一起。

第四，多切片可视化时，注意保持一致的坐标方向。不同切片在扫描时可能存在旋转或翻转，如果直接按原始像素坐标画图，切片间看起来方向不一致，影响对比。这种情况下，要么在实验阶段统一切片摆放方向，要么在可视化阶段对坐标做旋转或镜像变换。

提示：图像和坐标匹配还有一个容易忽略的点——如果实验流程对切片做了手工裁剪，那么图像中的组织轮廓和原始坐标文件可能就不再精确对应。这种情况下，优先相信坐标文件，以 spot 的相对位置关系为准，不要强求图像和 spot 像素级对齐。

最后再分享一个我自己的习惯。预处理和可视化跑完之后，我会把每个样本的关键 QC 指标汇总成一张表——总 spot 数、过滤后 spot 数、中位 UMI、中位基因数、线粒体比例中位数，再配上两到三张空间 QC 图。这张表和图组不仅是后续分析的基础，也是跟别人协作时最直观的交流材料。数据分析的结论可以被质疑，但完整的预处理记录和质控可视化是任何人都能查证的，这也是我做了这么多年数据分析最深的体会。