4. 数据加载与预处理：用 Scanpy/Seurat 搞定空间转录组数据
好，咱们正式开始动手干活了。数据加载和预处理，说白了就是给后续分析打好地基。我见过太多同学一上来就急着看聚类结果，结果因为预处理这步没做好，后面怎么调参数都怪怪的。今天我就把我在多个项目里踩过的坑和总结的经验，全盘托出。
4.1 环境准备：别在工具上浪费时间
我个人习惯用 conda 管理环境，干净、隔离、好复现。你想想看，要是因为包版本冲突把半天时间搭进去，多不值。
# 创建 scanpy 环境
conda create -n spa_scanpy python=3.10
conda activate spa_scanpy
pip install scanpy[leiden] anndata matplotlib pandas
# 或者创建 seurat 环境（R）
conda create -n spa_seurat r-base=4.3
conda activate spa_seurat
R -e "install.packages('Seurat')"
嗯，这里要注意：Scanpy 和 Seurat 的环境最好分开。我之前在同一个环境里混着装，结果 scanpy 的 matplotlib 和 seurat 的 ggplot2 互相干扰，排查了半天。
4.2 加载空间转录组数据：从原始格式到 AnnData/SeuratObject
空间转录组数据的格式五花八门。10x Visium 的数据最常见，包含三个核心文件：
- filtered_feature_bc_matrix.h5：表达矩阵
- spatial/tissue_positions_list.csv：spot 坐标
- spatial/scalefactors_json.json：缩放因子
我用 Scanpy 为例，给你演示一下标准加载流程。
import scanpy as sc
import pandas as pd
import numpy as np
# 设置全局参数
sc.settings.verbosity = 2 # 让日志详细点，方便排错
sc.settings.figdir = './figures/'
# 加载单样本 10x Visium 数据
adata = sc.read_visium(
path='./data/sample1/',
count_file='filtered_feature_bc_matrix.h5',
load_images=True
)
# 检查数据基本结构
print(adata)
print(adata.obs.head())
print(adata.obsm['spatial'][:5]) # 坐标信息
为什么会报错？最常见的原因是路径问题。我吃过一次亏：`count_file` 参数默认找 `raw_feature_bc_matrix.h5`，但空间转录组项目通常给的是过滤过的矩阵。所以记得指定文件名。
4.3 合并多个样本：不是简单的堆叠
多样本合并，看起来就是把几个 adata 拼接起来。但直接 `concat` 会丢失样本信息。我建议用两种方式：
| 方法 | 适用场景 | 关键参数 |
|---|---|---|
| 按批次合并（batch-aware） | 技术批次效应明显 | batch_key='sample_id' |
| 无批次合并 | 同一实验批次，处理一致 | index_unique='-' |
我个人习惯给每个样本先标记好批次信息，再合并。你看：
# 假设你有 3 个样本
samples = ['sample1', 'sample2', 'sample3']
adatas = []
for i, sample in enumerate(samples):
adata = sc.read_visium(f'./data/{sample}/')
adata.obs['sample_id'] = sample
adata.obs['batch'] = f'batch_{i}'
adatas.append(adata)
# 合并
adata_merged = sc.concat(adatas,
axis=0,
join='outer', # 保留所有基因
label='sample_id',
index_unique='-')
# 检查合并结果
print(f"合并后细胞数: {adata_merged.n_obs}")
print(f"合并后基因数: {adata_merged.n_vars}")
print(adata_merged.obs['sample_id'].value_counts())
注意：`join='outer'` 会保留所有样本中出现过的基因。如果你的数据是同一个平台生产的，基本没问题。但如果样本来自不同批次，里面有些低表达基因只在某个样本中出现，可能会带来噪声。
4.4 质量过滤：去掉那些不好看的 spot
合并后，第一件事就是过滤。你想想看，一个 spot 里只有几十个 UMI，或者基因数少得可怜，那很可能是个破损的组织边缘或者空spot。留着它，会让后续的归一化和聚类都跑偏。
我通常分两步走：
- 计算 QC 指标：UMI 总数、基因数、线粒体基因比例
- 设置过滤阈值：基于统计分布手动或自动设定
# 计算 QC 指标
adata_merged.var['mt'] = adata_merged.var_names.str.startswith('MT-')
sc.pp.calculate_qc_metrics(adata_merged, qc_vars=['mt'])
# 检查分布
sc.pl.violin(adata_merged,
['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],
groupby='sample_id',
multi_panel=True)
嗯，看小提琴图的时候，我一般关注两点：一是不同样本间的分布是否一致，二是尾部有没有明显离群点。如果一个样本的中位数远低于其他样本，说明那个样本质量可能有问题。
过滤代码也很直接：
# 手动设置阈值（根据提琴图调整）
sc.pp.filter_cells(adata_merged, min_genes=200)
sc.pp.filter_cells(adata_merged, max_genes=6000)
sc.pp.filter_genes(adata_merged, min_cells=3)
# 过滤高线粒体比例 spot
adata_merged = adata_merged[adata_merged.obs['pct_counts_mt'] < 20].copy()
# 过滤后看看还剩下多少
print(f"过滤后细胞数: {adata_merged.n_obs}")
print(f"过滤后基因数: {adata_merged.n_vars}")
4.5 基因筛选：保留有信息量的基因
过滤完 spot，下一步是挑基因。空间转录组通常有 2-5 万个基因，但大多数是低表达的 housekeeping 基因或者噪音。我们只保留高变基因，后续的降维和聚类才会聚焦于有生物学意义的变化。
Scanpy 里用 seurat_v3 方法效果不错，我实测过：
# 标准化对数化（seurat_v3 要求先做这个）
sc.pp.normalize_total(adata_merged, target_sum=1e4)
sc.pp.log1p(adata_merged)
# 挑选高变基因
sc.pp.highly_variable_genes(
adata_merged,
n_top_genes=2000,
batch_key='sample_id', # 按批次分别挑选
flavor='seurat_v3'
)
# 看看选了哪些
print(f"高变基因数: {adata_merged.var['highly_variable'].sum()}")
print(adata_merged.var[adata_merged.var['highly_variable']].head())
为什么用 seurat_v3 而不是默认的 seurat？因为前者考虑了方差与均值的关系，对低表达基因更友好。空间转录组的数据稀疏性比单细胞更高，低表达基因可能也很重要。
batch_key 可以避免某个样本的 batch effect 主导基因选择。如果不指定，一个样本里的高表达基因可能会挤掉其他样本的标记基因。
最后，只保留高变基因用于后续分析：
# 保存原始全部基因，只取高变基因继续
adata_merged.raw = adata_merged
adata_merged = adata_merged[:, adata_merged.var['highly_variable']].copy()
print(f"最终用于分析的数据: {adata_merged.n_obs} spots, {adata_merged.n_vars} genes")
4.6 本章知识结构：数据加载与预处理全景图
我把整个过程画成一张图，你看一眼就能记住关键节点。
你看，整个流程很清晰：从左到右是数据加载和合并，然后向下分流到过滤和筛选，最后汇聚成干净的数据集。每一步都有它的作用，跳一步都不行。