4. 数据读取与质量控制（QC）参数设置
做 Visium 分析，第一步就是把手里的原始数据喂给模型。这一步看着简单，但我见过太多人在这个环节翻车。数据没对齐、基因名格式不对、空点没过滤干净…这些问题到了聚类阶段全都会放大，最后结果根本没法看。
质量控制，说白了就是给数据“体检”。把那些明显有问题的点、噪音基因、背景污染剔除掉，咱们只保留高质量数据进行后续分析。这一步的参数设置，直接决定了聚类结果的下限。
4.1 核心数据读取函数与参数
我个人习惯用 sc.read_visium() 或 squidpy.read.visium() 来加载数据。这两个函数底层逻辑一样，都是解析 Space Ranger 输出的标准文件夹结构。
看一下典型用法：
import scanpy as sc
import squidpy as sq
# 读取原始数据，count_file 指定矩阵文件
adata = sc.read_visium(
path="./data/151673/", # Space Ranger 输出文件夹
count_file="filtered_feature_bc_matrix.h5",
load_images=True # 同时加载组织图像信息
)
# 或者用 Squidpy 版本
adata = sq.read.visium(
path="./data/151673/",
library_id="V10B01-084_A1", # 显式指定文库ID
source_image_path="./data/151673/spatial/tissue_hires_image.png"
)这里有几个关键点要注意：
- count_file 选项：Space Ranger 默认输出两个 h5 文件。
filtered_feature_bc_matrix.h5是已经过滤了空液滴的版本。raw_feature_bc_matrix.h5包含所有 barcode。我的建议：直接用 filtered 版本就行，少折腾。 - library_id 参数：如果你一个项目里有多张切片，一定要给每张切片指定唯一的 ID。不然后续拼接或对比分析时，数据会串。
- load_images 参数：设为
True才能拿到组织图像坐标。这个信息对于空间特征提取很关键。
adata 的概要信息，确认行数和列数是否符合预期。
4.2 数据文件的检验与处理
数据读进来以后，别急着往下跑。先做几项基本检验。我自己踩过坑：有一回跑聚类，结果死活不对，折腾了半天才发现是数据里混入了几千个空的 spot 坐标。
# 检验数据结构
print(adata)
print(adata.obs.head()) # 查看细胞/spot 信息
print(adata.var.head()) # 查看基因信息
print(adata.obsm.keys()) # 查看空间坐标等多元数据常见的数据格式问题和解决办法：
| 问题类型 | 表现 | 解决办法 |
|---|---|---|
| 基因名格式不标准 | 基因名为 ENSEMBL ID 而非符号 | 用 sc.queries.biomart_annotations() 转换 |
| 存在重复基因名 | 聚类特征高维冗余 | 用 adata.var_names_make_unique() 去重 |
| 空 spot 未过滤干净 | 大量零计数 barcode | 检查 adata.obs['total_counts'] 分布，手动过滤 |
4.3 质控参数设置详解
质量控制的核心思路就一条：识别并移除那些明显异常的数据点。我说“异常”，一般是指以下三种情况：
- 总基因数（UMI）过少的 spot —— 可能是空液滴或组织边缘。
- 总基因数过多的 spot —— 可能是两个 spot 融合，或区域组织破裂。
- 线粒体基因占比过高的 spot —— 细胞状态很差，可能已经死亡。
4.3.1 基因数（n_genes_by_counts）过滤
这个参数是最基础的。空间转录组的每个 spot 里理论上只有几个细胞。如果一个 spot 检测到的基因数特别少，它很可能就是空 spot 或背景区域。
# 计算每个 spot 的基因数
sc.pp.calculate_qc_metrics(adata, percent_top=None, log1p=False, inplace=True)
# 查看基因数分布
import matplotlib.pyplot as plt
plt.figure()
sc.pl.violin(adata, 'n_genes_by_counts', groupby='library_id')
plt.show()
# 设置过滤阈值，保留基因数 > 200 且 < 5000 的 spot
adata = adata[adata.obs['n_genes_by_counts'] < 5000, :].copy()
adata = adata[adata.obs['n_genes_by_counts'] > 200, :].copy()阈值怎么定？没有标准答案。我的经验是：先画个小提琴图看看分布，下界选在“明显掉头”的位置，上界去掉最右边 2% 的离群值。不同组织差异很大，肿瘤样本的 spot 基因数普遍偏高。
4.3.2 总UMI数（total_counts）与线粒体比例
基因数考虑的是“有多少种基因被检测到了”，总 UMI 数考虑的是“总共检测到了多少条转录本”。两者侧重点不一样。线粒体比例则是“红线指标”，一旦偏高，这个 spot 的数据基本不可信。
# 计算线粒体基因比例
adata.var['mt'] = adata.var_names.str.startswith('MT-') # 人源
# 或 adata.var['mt'] = adata.var_names.str.startswith('mt-') # 鼠源
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)
# 过滤：总 UMI 数 500~30000，线粒体比例 < 20%
adata = adata[adata.obs['total_counts'] > 500, :].copy()
adata = adata[adata.obs['total_counts'] < 30000, :].copy()
adata = adata[adata.obs['pct_counts_mt'] < 20, :].copy()线粒体比例的上限，我一般设在 20% 到 25%。如果组织样本不太新鲜，比例可能会偏高，可以放宽到 30%。但超过 30% 的 spot 我建议直接丢弃，留着只会拉低聚类信噪比。
4.3.3 细胞周期评分（可选）
有些场景下，我会额外加一个细胞周期评分的质控步骤。这步不是必须的，但对于研究增殖活跃的组织（比如肿瘤或发育组织），它能帮你识别出“假阳性高表达”的 spot。
# 使用 Scanpy 内置的细胞周期基因集（人源）
cell_cycle_genes = [gene for gene in adata.var_names if 'S' in gene or 'G2M' in gene]
sc.tl.score_genes_cell_cycle(adata, s_genes=s_genes, g2m_genes=g2m_genes)
# 查看细胞周期评分的分布，异常高值可能提示分裂中的细胞团
sc.pl.violin(adata, ['S_score', 'G2M_score'])4.4 完整的质量控制流程
把上面的步骤串起来，就是一个完整的 QC 流水线。我建议把流程写成一个函数，方便复用在多张切片上。
def run_visium_qc(adata,
min_genes=200,
max_genes=5000,
min_umis=500,
max_umis=30000,
max_mito=0.20,
species='human'):
# Step 1: 基础统计
if species == 'human':
adata.var['mt'] = adata.var_names.str.startswith('MT-')
else:
adata.var['mt'] = adata.var_names.str.startswith('mt-')
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)
# Step 2: 按阈值过滤
mask = (adata.obs['n_genes_by_counts'] > min_genes) & \
(adata.obs['n_genes_by_counts'] < max_genes) & \
(adata.obs['total_counts'] > min_umis) & \
(adata.obs['total_counts'] < max_umis) & \
(adata.obs['pct_counts_mt'] < max_mito * 100)
adata_filtered = adata[mask].copy()
# Step 3: 标准化（为了下游分析）
sc.pp.normalize_total(adata_filtered, target_sum=1e4)
sc.pp.log1p(adata_filtered)
print(f"QC 保留 spots: {adata_filtered.shape[0]} / {adata.shape[0]}")
return adata_filtered调用起来很简单：
adata_qc = run_visium_qc(adata, min_genes=100, max_genes=4000, max_mito=0.15)4.5 我应该如何调整这些参数？
这可能是大家最关心的问题。我总结了几种典型场景下的参考配置，但记住：每个数据集都不一样，一定要可视化确认后再批量跑。
| 数据类型 | min_genes | max_genes | max_mito | 说明 |
|---|---|---|---|---|
| 高质量新鲜组织 | 200 | 6000 | 0.20 | 最常用，平衡性最好 |
| FFPE 组织 | 100 | 3000 | 0.30 | RNA降解严重，适当放宽线粒体 |
| 肿瘤样本 | 300 | 10000 | 0.15 | 高异质性+高表达，收紧上界再降低线粒体 |
| 发育组织（胚胎） | 150 | 5000 | 0.25 | 细胞类型复杂，基因数差异大 |
另外，我强烈建议你画一张 QC 指标的联合分布图。它可以直观显示过滤前后数据的变化：
sc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts', color='pct_counts_mt')从图上你一眼就能看出异常点在哪。那些“矮胖”的点往往组织质量差，早点去掉免得干扰聚类。
4.6 避坑指南：我经历过的三个 QC 惨案
嗯，既然聊到参数设置，我就说说我犯过的错。希望你能绕开这些坑。
- 惨案一：忽略组织图像对齐 — 有一次我直接跑默认读取，结果数据矩阵和图像坐标错位了。聚类结果一路飘红，最后才发现是初始化路劲里
library_id没对应上。解决方案：读取后一定要验证adata.obsm['spatial']的坐标是否和 H&E 图像对齐。 - 惨案二：线粒体阈值设太死 — 在 FFPE 样本上用 15% 的线粒体上限，结果删掉了 70% 的 spot。你想想看，一个切片只剩 200 多个 spot，后面还怎么聚类？所以对不同样本类型，阈值一定要灵活。
- 惨案三：盲目套用单细胞 RNA-seq 的 QC 标准 — 空间转录组和单细胞不一样，每个 spot 有多个细胞混合，基因数的分布范围宽得多。不要用 scRNA-seq 那套严格的 200~2500 基因数限制，否则大多数 spot 都会被过滤掉。
4.7 小结：质控参数的平衡艺术
质量控制不是越严格越好，也不是越宽松越好。关键是在保持数据量的同时，尽可能剔除噪音。我每次做一个新数据集，都会先跑 3~4 组不同的 QC 参数做对比，再决定最终方案。
你可以把今天学到的流程固化下来，做成一个可复用的脚本。未来每拿到一张新切片，直接跑一遍，看看需要微调哪些参数。相信我，这会节省你大量调试时间。
_raw 和 _qc 后缀。