一、前言

空间转录组技术（Spatial Transcriptomics）可以同时获取组织切片中每个空间位置的基因表达信息，具有重要的研究价值，如用于组织结构分区、细胞类型识别、病理状态分析等。
本教程基于 Scanpy，结合 10X Genomics 提供的空间转录组数据格式，展示完整的数据读取、预处理、质控、可视化、聚类及差异分析流程。

二 、代码部分

1. 读取数据

构建转录组数据的目录结构

library_id = 'BC_06066'
data_dir = r'F:\algorithm\python\feature-mapping\mapping_related'
ST_dir = os.path.join(data_dir, 'ST')
BC_dir = os.path.join(ST_dir, library_id)

h5_path = os.path.join(BC_dir, 'filtered_feature_bc_matrix.h5')
spatial_dir = os.path.join(BC_dir, 'spatial')
position_path = os.path.join(spatial_dir, 'tissue_positions_list.csv')
hires_image_path = os.path.join(spatial_dir, 'tissue_hires_image.png')
scalefactors_json_path = os.path.join(spatial_dir, 'scalefactors_json.json')

这些路径构建了标准的空间转录组数据结构，根据自己数据的实际情况进行构建。我的数据包括：

h5_path：表达矩阵路径（.h5 文件）；
position_path：空间点位信息；
hires_image_path：原始高分辨率病理组织图像；
scalefactors_json_path：缩放因子配置。

读取表达矩阵并处理变量名

adata = sc.read_10x_h5(h5_path)
adata.var_names_make_unique()
print(adata)

2. 加载空间坐标

positions = pd.read_csv(position_path, header=None)
positions.columns = ["barcode", "in_tissue", "array_row", "array_col", "pxl_row_in_fullres", "pxl_col_in_fullres"]
positions.set_index('barcode', inplace=True)

adata.obs = adata.obs.join(positions, how='left')
adata.obsm['spatial'] = adata.obs[['pxl_row_in_fullres', 'pxl_col_in_fullres']].to_numpy()

加载组织图像与缩放因子

hires_image = plt.imread(hires_image_path)
with open(scalefactors_json_path, 'r') as f:
    scalefactors = json.load(f)

adata.uns['spatial'] = {
    library_id: {
        'images': {
            'hires': hires_image
        },
        'scalefactors': scalefactors
    }
}

3. 数据预处理与质控

基本过滤

sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)

过滤掉基因数低于 200 的细胞；
去除在少于 3 个细胞中表达的基因。

添加线粒体基因标识并计算 QC 指标

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)

添加布尔型 mt 字段，标记线粒体基因；
计算细胞总转录本数、基因数、线粒体基因占比等质控指标。

筛选高质量细胞

adata = adata[adata.obs.pct_counts_mt < 5, :]
adata = adata[adata.obs['in_tissue'] == 1, :]

保留线粒体基因比例低于 5% 的细胞；
保留在组织切片区域内的细胞（in_tissue == 1）。

归一化与对数变换

sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

对每个细胞的表达量总和进行归一化（总和为 10,000）；
使用自然对数进行平滑，降低极端表达值的影响。

4. 空间表达可视化

sc.pl.spatial(adata, img_key="hires", color="total_counts", frameon=False, title="BC_06066 - Total UMI", cmap="viridis")
sc.pl.spatial(adata, img_key="hires", color="n_genes", frameon=False, title="BC_06066 - Total Gene", cmap="Spectral_r")

展示细胞在组织切片上的分布情况，配色反映 UMI 数或基因数量。

5. 高变基因筛选与标准化处理

筛选高变基因

sc.pp.highly_variable_genes(adata, flavor='seurat', n_top_genes=2000)
adata = adata[:, adata.var['highly_variable']]

使用 Seurat 方法选出 2000 个变异性最大的基因用于后续分析，有助于聚类与降维。

回归处理与标准化

sc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt'])
sc.pp.scale(adata, max_value=10)

回归掉测序深度和线粒体占比；
将每个基因缩放为均值为 0、标准差为 1，并截断至 [-10, 10] 范围。

6. PCA 降维与聚类

PCA 主成分分析

sc.tl.pca(adata, n_comps=50, svd_solver='arpack')
sc.pl.pca_variance_ratio(adata, log=True, n_pcs=50)

使用 PCA 将表达矩阵压缩为 50 个主成分；
可视化主成分的解释度，评估降维效果。

构建邻居图并leiden聚类

sc.pp.neighbors(adata, n_neighbors=20, n_pcs=50)
sc.tl.leiden(adata)
sc.pl.spatial(adata, img_key="hires", color="leiden", frameon=False)

基于 PCA 结果构建 KNN 图；
使用 Leiden 算法聚类，并将聚类结果在空间中可视化。

7. 非线性降维：t-SNE 和 UMAP

sc.tl.tsne(adata)
sc.pl.tsne(adata, color='leiden')

sc.tl.umap(adata)
sc.pl.umap(adata, color='leiden')

t-SNE 和 UMAP 是常用的非线性降维方法，适用于高维转录组数据的可视化展示聚类分布。

8. 差异基因分析

sc.tl.rank_genes_groups(adata, 'leiden', method='t-test')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)

result = adata.uns['rank_genes_groups']
groups = result['names'].dtype.names

for group in groups:
    print(f"Cluster {group}:")
    for gene in result['names'][group][:5]:
        print(f"  {gene}")

对每个聚类（leiden）执行差异基因分析，找出代表性 marker；
默认方法为 t-test，可切换为 wilcoxon 或 logreg。

9. 聚类标签命名与可视化

new_cluster_names = ['IGHG2', 'PPAN', 'SPR', 'APOD', 'BGN', 'C1R', 'C3', 'CAV1', 'EGR1', 'FAS', 'LRP1', 'PER1']
adata.rename_categories('leiden', new_cluster_names)
sc.pl.umap(adata, color='leiden', legend_loc='on data', title='', frameon=False, save='.pdf')

将 Leiden 聚类结果替换为具有生物学意义的基因名，便于进一步分析解释。

10. 完整代码

import json
import scanpy as sc
import matplotlib.pyplot as plt
import pandas as pd
import os
'''
读取数据
'''
library_id = 'BC_06066'
data_dir = r'F:\algorithm\python\feature-mapping\mapping_related'  # 根目录
ST_dir = os.path.join(data_dir, 'ST')  # ST目录

BC_dir = os.path.join(ST_dir, library_id)  # ST其中1例ST数据
h5_path = os.path.join(BC_dir, 'filtered_feature_bc_matrix.h5')  # h5文件路径
spatial_dir = os.path.join(BC_dir, 'spatial')  # 空间信息路径
position_path = os.path.join(spatial_dir, 'tissue_positions_list.csv')  # 空间坐标信息路径
hires_image_path = os.path.join(spatial_dir, 'tissue_hires_image.png')  # 高分辨率图像路径
scalefactors_json_path = os.path.join(spatial_dir, 'scalefactors_json.json')  # 缩放因子路径

# 读取基因表达矩阵，读取未过滤的自己进行质控。
adata = sc.read_10x_h5(h5_path)
adata.var_names_make_unique()
print(adata)

'''
obs：细胞相关的元数据（行索引是细胞ID）。
var：基因相关的元数据（行索引是基因ID）。
uns：非结构化的注释数据，存储各种辅助信息和参数。
obsm：细胞级别的多维数据，例如空间坐标。
'''

# 读取空间坐标表
positions = pd.read_csv(position_path, header=None)
positions.columns = ["barcode", "in_tissue", "array_row", "array_col", "pxl_row_in_fullres", "pxl_col_in_fullres"]
# 将barcode设置为索引
positions.set_index('barcode', inplace=True)
# 将空间坐标信息左连接添加到adata.obs中
adata.obs = adata.obs.join(positions, how='left')
# 将空间坐标添加到adata.obsm中
adata.obsm['spatial'] = adata.obs[['pxl_row_in_fullres', 'pxl_col_in_fullres']].to_numpy()

# 其他的空间信息在adata.uns中
hires_image = plt.imread(hires_image_path)  # 高分辨率图像地址
with open(scalefactors_json_path, 'r') as f:  # 加载json文件中的缩放因子
    scalefactors = json.load(f)
adata.uns['spatial'] = {
    library_id: {
        'images': {
            'hires': hires_image
        },
        'scalefactors': scalefactors
    }
}

'''
数据预处理与质控
'''
# 过滤掉那些表达的基因数量少于200个的低质量的细胞，添加细胞的表达基因数n_genes到obs中
sc.pp.filter_cells(adata, min_genes=200)
# 过滤掉那些在少于3个细胞中表达的基因，添加基因对应的细胞数n_cells到var中
sc.pp.filter_genes(adata, min_cells=3)
# 识别所有mt开头的基因即线粒体基因，添加mt属性(bool值)到var中，但给的数据里没有线粒体基因
adata.var['mt'] = adata.var_names.str.startswith('MT-')
# 计算质量控制（QC）指标，将结果属性添加到obs和var中
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)

# # 线粒体基因占比的可视化 - 小提琴图
# sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],jitter=0.4, multi_panel=True)
# # 线粒体基因占比的可视化 - 散点图
# sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt')
# sc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts')

# 过滤掉线粒体基因比例超过5%的细胞
adata = adata[adata.obs.pct_counts_mt < 5, :]
# 过滤掉组织外的数据
adata = adata[adata.obs['in_tissue'] == 1, :]
# 进行数据的总量归一化，将每个细胞的总UMI数归一化为10,000。
sc.pp.normalize_total(adata, target_sum=1e4)
# 对数据进行对数变换，使得表达数据更接近正态分布
sc.pp.log1p(adata)

# # 绘制前20个表达最高的基因
# sc.pl.highest_expr_genes(adata, n_top=20, )

'''
空间位置可视化
'''
sc.pl.spatial(adata, img_key="hires", color="total_counts", frameon=False,
              title="BC_06066 - Total UMI", cmap="viridis")
sc.pl.spatial(adata, img_key="hires", color="n_genes", frameon=False,
              title="BC_06066 - Total Gene", cmap="Spectral_r")

'''
降维、聚类
'''
# 获取前2000个高变基因
sc.pp.highly_variable_genes(adata, flavor='seurat', n_top_genes=2000)
# 过滤非高变基因(可选)
adata = adata[:, adata.var['highly_variable']]

# 获取前20个高变基因打印名称
# top_var_genes = adata.var_names[adata.var.highly_variable.values][:20]
# print(top_var_genes)
# # 高变基因可视化
# sc.pl.highly_variable_genes(adata)

# 对数据进行回归处理，移除总UMI数量 (total_counts) 和线粒体基因比例 (pct_counts_mt) 对基因表达的影响。
sc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt'])
# 对每个基因的表达量进行标准化（均值为0，标准差为1），保留表达量在[-10, 10] 之间的数据。
sc.pp.scale(adata, max_value=10)

# PCA降维
sc.tl.pca(adata, n_comps=50, svd_solver='arpack')
sc.pl.pca_variance_ratio(adata, log=True, n_pcs=50) # 绘制前20个主成分的方差贡献率

# 图聚类算法leiden,并保存聚类结果到obs
sc.pp.neighbors(adata, n_neighbors=20, n_pcs=50)  # 根据pca结果构建邻居图
sc.tl.leiden(adata)
sc.pl.spatial(adata, img_key="hires", color="leiden", frameon=False)

# t-SNE降维
sc.tl.tsne(adata)
sc.pl.tsne(adata, color='leiden')

# UMAP降维
sc.tl.umap(adata)
sc.pl.umap(adata, color='leiden')

'''
差异性分析
'''
# 采用三种方法寻找25个maker基因分别是t-test、wilcoxon、logreg
sc.settings.verbosity = 2  # 设置日志输出

sc.tl.rank_genes_groups(adata, 'leiden', method='t-test')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)

# sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
# sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)
#
# sc.tl.rank_genes_groups(adata, 'leiden', method='logreg')
# sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)

# 查看排名前几的差异基因
result = adata.uns['rank_genes_groups']
groups = result['names'].dtype.names

# 打印每个聚类的前几个差异基因
for group in groups:
    print(f"Cluster {group}:")
    for gene in result['names'][group][:5]:  # 打印前5个基因
        print(f"  {gene}")

'''
确定相应的细胞类型，并分配给簇
'''

new_cluster_names = ['IGHG2', 'PPAN', 'SPR', 'APOD', 'BGN', 'C1R', 'C3', 'CAV1', 'EGR1',
                     'FAS', 'LRP1', 'PER1']
adata.rename_categories('leiden', new_cluster_names)
sc.pl.umap(adata, color='leiden', legend_loc='on data', title='', frameon=False, save='.pdf')

总结

本文通过 Scanpy 工具链，完成了从空间转录组数据读取、质控、归一化、降维、聚类、差异分析到细胞类型注释的完整流程，适用于基于 10X Visium 平台的空间数据分析场景。