单细胞空间转录组：基础数据分析（基于 Scanpy 的单样本流程）

作者: SeekGene

时长: 20 分钟

字数: 4.7k 字

更新: 2026-06-25

阅读: 0 次

空间转录组 Notebooks

1. 模块简介

SeekSpace 技术不仅能够检测到单细胞精度的基因表达数据，还能精准定位每个细胞在组织中的真实空间坐标，是一种高分辨率空间转录组技术。为了方便客户快速上手，SeekSpace 产生的数据能够无缝兼容主流的单细胞与空间转录组分析生态，例如 Python 的 scanpy 库和 R 的 Seurat 库。

本教程专为使用 Python 分析环境的客户编写。在这里，我们将演示如何生成和读取由 seekspacetools 生成的 .h5ad 文件，并进行直观的可视化分析。

2. 输入文件准备

原始输入数据说明： seekspacetools 也会生成矩阵结果文件目录，主要包含：

表达矩阵目录（如 *_filtered_feature_bc_matrix）：包含 barcodes.tsv.gz, features.tsv.gz, matrix.mtx.gz，记录了基因在每个细胞中的表达量；以及 cell_locations.tsv.gz，记录了细胞的物理空间坐标。
组织切片图像（如 *aligned_DAPI.png 和 *aligned_HE.png）：DAPI 为细胞核荧光染色图片，HE 为苏木精-伊红染色图片，用于展示组织形态学结构。

关于 .h5ad 文件： 这是一种业内标准的单细胞/空间转录组数据存储格式。官方的预处理流程已经为您完成了以下繁琐的早期工作：

数据的质控 (QC) 与细胞/基因过滤
表达量的归一化与高变基因 (HVG) 提取
降维与聚类 (PCA / UMAP / tSNE / Leiden 无监督分群)
空间坐标 (X, Y) 及组织染色图片 (DAPI / HE) 的提取与整合

这意味着，您只需要简单地加载该文件，即可直接进入高级的下游分析和可视化展示阶段！

💡 提示： SeekSpace 芯片中，一个像素点的大小约为 0.2653 微米。如果您需要计算细胞在真实组织空间上的物理距离，只需将像素坐标乘以 0.2653 即可得出实际的微米级距离。

python

import scanpy as sc
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import warnings
warnings.filterwarnings("ignore")

# 设置 scanpy 的日志显示级别和绘图全局参数
sc.settings.verbosity = 3
sc.settings.set_figure_params(dpi=100, facecolor="white")

3. 预处理流程：从原始矩阵到 h5ad

为了让您完全掌握数据的处理细节，本节展示了 seekspacetools 处理生成.h5ad文件的核心代码逻辑。它展示了如何从 seekspacetools 生成的表达矩阵、空间坐标文件以及组织图像，一步步构建出包含降维和聚类信息的完整 .h5ad 对象。

注意： 本节主要用于教学展示，帮助您理解底层逻辑。如果您只想快速进行可视化探索，可以跳过本节的代码，直接从第 4 节开始运行。

3.1 读取表达矩阵与空间坐标

python

# 1. 读取 CellRanger 格式的表达矩阵
adata_raw = sc.read_10x_mtx('./WTH1092_filtered_feature_bc_matrix/', var_names='gene_symbols', cache=False)
adata_raw.obs_names_make_unique()
adata_raw.var_names_make_unique()
# 
# 2. 读取细胞空间坐标文件 (cell_locations.tsv.gz)
df_center = pd.read_csv('./WTH1092_filtered_feature_bc_matrix/cell_locations.tsv.gz', sep='\t')
# 
# 3. 将空间坐标与表达矩阵的细胞对齐
df_center_aligned = df_center.set_index('Cell_Barcode').reindex(adata_raw.obs_names)
# 
# 4. Y轴坐标翻转 (根据图像真实物理尺寸 size_y 进行翻转，以适配 Scanpy 空间绘图系统)
size_y = 19906 # 示例图像Y轴真实尺寸
df_center_aligned["Y"] = size_y - df_center_aligned["Y"]
# 
# 5. 整合坐标到 adata 中，并保存为兼容 Seurat 习惯的 Spatial1/Spatial2 列
adata_raw.obsm['spatial'] = df_center_aligned[['X', 'Y']].to_numpy()
adata_raw.obs['Spatial1'] = df_center_aligned['X']
adata_raw.obs['Spatial2'] = df_center_aligned['Y']

output

--> This might be very slow. Consider passing \`cache=True\`, which enables much faster reading from a cache file.

3.2 质控 (QC) 与过滤

python

# 1. 过滤极低表达的细胞 (可选)
sc.pp.filter_genes(adata_raw, min_cells=0)
# 
# 2. 标记线粒体基因
adata_raw.var['mt'] = adata_raw.var_names.str.startswith(('MT-', 'mt-', 'Mt-'))
# 
# 3. 计算质控指标
sc.pp.calculate_qc_metrics(adata_raw, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)
# 
# 4. 重命名 QC 列名以兼容后续流程与 Seurat 习惯
adata_raw.obs.rename(columns={
    'n_genes_by_counts': 'nFeature_RNA',
    'total_counts': 'nCount_RNA',
    'pct_counts_mt': 'percent.mito'
}, inplace=True)
# 
# 5. 绘制质控小提琴图和散点图
sc.pl.violin(adata_raw, ['nFeature_RNA', 'nCount_RNA', 'percent.mito'], jitter=0.4, multi_panel=True, show=False)
sc.pl.scatter(adata_raw, x='nCount_RNA', y='percent.mito', show=False)
sc.pl.scatter(adata_raw, x='nCount_RNA', y='nFeature_RNA', show=False)

output

3.3 归一化、寻找高变基因 (HVG) 与数据标准化

python

# 1. 表达量归一化 (Normalize)
sc.pp.normalize_total(adata_raw, target_sum=1e4)
# 
# 2. 对数转换 (Log-transform)
sc.pp.log1p(adata_raw)
# 
# 3. 寻找前 2000 个高变基因 (使用 Seurat 算法)
sc.pp.highly_variable_genes(adata_raw, n_top_genes=2000, flavor='seurat')
# 
# 4. 冻结原始表达量以备后续差异分析使用
adata_raw.raw = adata_raw
# 
# 5. 提取高变基因子集，回归掉线粒体比例的影响，并进行标准化 (Scale)
adata_hvg = adata_raw[:, adata_raw.var.highly_variable].copy()
sc.pp.regress_out(adata_hvg, ['percent.mito'])
sc.pp.scale(adata_hvg, max_value=10)

output

normalizing counts per cell
finished (0:00:00)
extracting highly variable genes
finished (0:00:00)
--> added
'highly_variable', boolean vector (adata.var)
'means', float vector (adata.var)
'dispersions', float vector (adata.var)
'dispersions_norm', float vector (adata.var)
regressing out ['percent.mito']
sparse input is densified and may lead to high memory use
finished (0:01:18)

3.4 降维、无监督聚类与差异基因分析

python

# 1. 主成分分析 (PCA)
sc.tl.pca(adata_hvg, svd_solver='arpack')
adata_raw.obsm['X_pca'] = adata_hvg.obsm['X_pca']
# 
# 2. 构建细胞近邻图 (Neighborhood graph)
sc.pp.neighbors(adata_raw, n_neighbors=10, n_pcs=15)
# 
# 3. 遍历多个分辨率进行 Leiden 聚类 (0.2 ~ 1.7)
import numpy as np
for res in np.arange(0.2, 1.7, 0.3):
    sc.tl.leiden(adata_raw, resolution=res, key_added=f'res_{res:.1f}')
# 
# 4. 设定默认聚类结果 (resolution=0.8)
sc.tl.leiden(adata_raw, resolution=0.8, key_added='leiden')
# 
# 5. 非线性降维 (UMAP 与 tSNE)
sc.tl.umap(adata_raw)
sc.tl.tsne(adata_raw, n_pcs=15)
# 
# 6. 寻找差异表达基因 (Marker genes)
sc.tl.rank_genes_groups(adata_raw, 'leiden', method='wilcoxon', key_added='rank_genes_groups')

output

computing PCA
with n_comps=50
finished (0:00:12)
computing neighbors
using 'X_pca' with n_pcs = 15
finished: added to \`.uns['neighbors']\`
\`.obsp['distances']\`, distances for each pair of neighbors
\`.obsp['connectivities']\`, weighted adjacency matrix (0:00:04)
running Leiden clustering
finished: found 12 clusters and added
'res_0.2', the cluster labels (adata.obs, categorical) (0:00:20)
running Leiden clustering
finished: found 14 clusters and added
'res_0.5', the cluster labels (adata.obs, categorical) (0:00:18)
running Leiden clustering
finished: found 18 clusters and added
'res_0.8', the cluster labels (adata.obs, categorical) (0:00:24)
running Leiden clustering
finished: found 23 clusters and added
'res_1.1', the cluster labels (adata.obs, categorical) (0:00:19)
running Leiden clustering
finished: found 27 clusters and added
'res_1.4', the cluster labels (adata.obs, categorical) (0:00:27)
running Leiden clustering
finished: found 18 clusters and added
'leiden', the cluster labels (adata.obs, categorical) (0:00:24)
computing UMAP
finished: added
'X_umap', UMAP coordinates (adata.obsm)
'umap', UMAP parameters (adata.uns) (0:00:19)
computing tSNE
using 'X_pca' with n_pcs = 15
using sklearn.manifold.TSNE
finished: added
'X_tsne', tSNE coordinates (adata.obsm)
'tsne', tSNE parameters (adata.uns) (0:02:28)
ranking genes
finished: added to \`.uns['rank_genes_groups']\`
'names', sorted np.recarray to be indexed by group ids
'scores', sorted np.recarray to be indexed by group ids
'logfoldchanges', sorted np.recarray to be indexed by group ids
'pvals', sorted np.recarray to be indexed by group ids
'pvals_adj', sorted np.recarray to be indexed by group ids (0:01:38)

3.5 整合组织切片图像并保存 h5ad

python

import os
from PIL import Image
# 
# 1. 读取组织切片图片 (DAPI和HE) 并根据图像像素宽度计算真实的 scalefactor
img_dapi_path = "WTH1092_aligned_DAPI.png"
img_he_path = "WTH1092_aligned_HE.png"

# size_y和size_y保存在 h5ad文件的 adata.uns['spatial'][sample_name] 中。
## 如您的软件版本为 v1.0.2, size_x 和 size_y 如下
size_x = 55050 
size_y = 19906 
## 如您的软件版本为 v1.0.3, size_x 和 size_y 如下
# size_x = 56500 
# size_y = 20434 


tissue_DAPI_scalef = 1.0
tissue_HE_scalef = 1.0
# 
img_dapi = Image.open(img_dapi_path).convert("RGB")
img_dapi_width, _ = img_dapi.size
tissue_DAPI_scalef = img_dapi_width / size_x
# 
if os.path.exists(img_he_path):
    img_he = Image.open(img_he_path).convert("RGB")
    img_he_width, _ = img_he.size
    tissue_HE_scalef = img_he_width / size_x
else:
    tissue_HE_scalef = tissue_DAPI_scalef
# 
# 2. 将图像矩阵及空间元数据存入 adata.uns['spatial'] 中
library_id = "WTH1092"
adata_raw.uns['spatial'] = {
    library_id: {
        "images": {"DAPI": np.array(img_dapi)},
        "scalefactors": {
            "tissue_DAPI_scalef": tissue_DAPI_scalef,
            "tissue_HE_scalef": tissue_HE_scalef,
            "spot_diameter_fullres": 100.0
        },
        "size_x": size_x,
        "size_y": size_y
    }
}
if os.path.exists(img_he_path):
    adata_raw.uns['spatial'][library_id]['images']['HE'] = np.array(img_he)
# 
# 3. 将处理好的完整数据对象保存为压缩的 h5ad 文件，供下游分析直接加载！
adata_raw.write_h5ad("WTH1092.h5ad", compression="gzip")

4. 读取 h5ad 综合数据对象 (快速开始)

请将下方代码中的 'B21_36M.h5ad' 替换为您实际生成的 .h5ad 样本文件路径。读取后，您将得到一个 AnnData 对象（通常命名为 adata），它是贯穿整个 Python 空间转录组分析的核心数据结构。

python

# 读取预处理后的综合数据文件
# 这里的 'B21_36M.h5ad' 请根据您的实际文件名进行修改
adata = sc.read_h5ad('WTH1092.h5ad')

# 查看 adata 对象的概览信息
# 您可以在输出中看到细胞数(n_obs)、基因数(n_vars) 以及存储的降维和空间信息
adata

output

AnnData object with n_obs × n_vars = 32758 × 32285
obs: 'Spatial1', 'Spatial2', 'nFeature_RNA', 'nCount_RNA', 'total_counts_mt', 'percent.mito', 'res_0.2', 'res_0.5', 'res_0.8', 'res_1.1', 'res_1.4', 'leiden'
var: 'gene_ids', 'feature_types', 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'highly_variable', 'means', 'dispersions', 'dispersions_norm'
uns: 'hvg', 'leiden', 'log1p', 'neighbors', 'rank_genes_groups', 'res_0.2', 'res_0.5', 'res_0.8', 'res_1.1', 'res_1.4', 'spatial', 'tsne', 'umap'
obsm: 'X_pca', 'X_tsne', 'X_umap', 'spatial'
obsp: 'connectivities', 'distances'

5. UMAP 降维与细胞分群展示

在转录组分析中，我们通常会将高维的基因表达数据降维到二维平面（如 UMAP 或 tSNE）进行展示。在降维图上，距离越近的细胞，其整体基因表达模式越相似。

官方预处理脚本已经利用无监督聚类算法（Leiden）将细胞分为了不同的亚群（Cluster），并将结果存储在 adata.obs['leiden'] 中。

python

# 设置绘图尺寸
sc.settings.set_figure_params(figsize=(7, 7))

# 按照 Leiden 聚类结果在 UMAP 图上进行着色
sc.pl.umap(adata, color=['leiden'])

除了展示细胞群体，您还可以查看特定感兴趣的 Marker 基因在各个细胞群中的表达丰度分布：

python

# 这里以 'Cst3' 基因作为示例，您可以将其替换为您研究的关键基因名称
gene_of_interest = 'Cst3'
if gene_of_interest in adata.var_names:
    sc.pl.umap(adata, color=[gene_of_interest], use_raw=False, cmap='Reds')
else:
    print(f"未在数据中找到基因 {gene_of_interest}")

6. 原位空间分布展示

SeekSpace 的核心技术优势在于保留了细胞在组织切片上的真实物理位置！

细胞的空间坐标信息已被安全地存储在 adata.obsm['spatial'] 中。我们可以通过 scanpy 的 embedding 可视化功能，直接展示细胞在芯片空间上的分布，并利用不同的颜色来标记它们所属的细胞类群。

python

sc.settings.set_figure_params(figsize=(12, 5))

# 绘制原位空间散点图，颜色映射为细胞所属的分群
sc.pl.embedding(adata, "spatial", color=["leiden"])

7. 叠加组织切片图像展示

为了更加直观地将细胞类群与真实的组织病理学结构对应，我们可以将细胞的空间散点图直接覆盖在组织的荧光染色（DAPI）或明场染色（HE）图片上。

官方脚本已经为您提取并保存了这些背景图像及精确的缩放比例（Scalefactors），您可以使用 sc.pl.spatial 函数一键绘制出兼具科研严谨性与美观度的空间映射图。

python

# 获取当前数据的空间库名称 (library_id)
library_id = list(adata.uns['spatial'].keys())[0]

# 优先尝试使用 DAPI 图像作为背景
if 'DAPI' in adata.uns['spatial'][library_id]['images']:
    sc.pl.spatial(adata, color="leiden", library_id=library_id, img_key="DAPI", 
                  size=1.5, alpha=1.0, title="Spatial distribution overlaid on DAPI")

# 如果样本提供了 HE 图像，也可以使用 HE 图像作为背景
if 'HE' in adata.uns['spatial'][library_id]['images']:
    sc.pl.spatial(adata, color="leiden", library_id=library_id, img_key="HE", 
                  size=1.5, alpha=1.0, title="Spatial distribution overlaid on HE")

7.1 (进阶) 使用 Matplotlib 高度自定义空间图

如果您在准备学术论文插图或汇报时，需要对图像的细节（如细胞点的大小、颜色透明度、图例位置、坐标轴范围等）进行更深度的个性化调整，可以参考以下使用 matplotlib 手动叠加背景图和坐标点的代码模板：

python

# 提取组织图像和真实的尺寸范围信息
if 'DAPI' in adata.uns['spatial'][library_id]['images']:
    background_image = adata.uns['spatial'][library_id]['images']['DAPI']
elif 'HE' in adata.uns['spatial'][library_id]['images']:
    background_image = adata.uns['spatial'][library_id]['images']['HE']
else:
    background_image = None

size_x = adata.uns['spatial'][library_id]['size_x']
size_y = adata.uns['spatial'][library_id]['size_y']

# 提取空间坐标和细胞分群信息到一个方便操作的 DataFrame 中
spatial_df = pd.DataFrame(adata.obsm['spatial'], columns=['spatial_1', 'spatial_2'], index=adata.obs.index)
spatial_df['leiden'] = adata.obs['leiden']

# 预设一些科研绘图中常用且区分度高的美观配色
MYCOLOR=[
    "#6394ce", "#2a4c87", "#eed500", "#ed5858", "#f6cbc2", "#f5a2a2", "#3ca676", "#6cc9d8",
    "#ef4db0", "#992269", "#bcb34a", "#74acf3", "#3e275b", "#fbec7e", "#ec4d3d", "#ee807e",
    "#f7bdb5", "#dbdde6", "#f591e1", "#51678c", "#2fbcd3", "#80cfc3", "#fbefd1", "#edb8b5",
    "#5678a8", "#2fb290", "#a6b5cd", "#90d1c1", "#a4e0ea", "#837fd3", "#5dce8b", "#c5cdd9",
    "#f9e2d6", "#c64ea4", "#b2dfd6", "#dbdfe7", "#dff2ec", "#cce8f3", "#e74d51", "#f7c9c4"
]
cluster = list(set(adata.obs["leiden"]))
colors = {cluster[i]: MYCOLOR[i % len(MYCOLOR)] for i in range(len(cluster))}

# ---------------- 开始绘制 ----------------
fig, ax = plt.subplots(figsize=(12, 8))
if background_image is not None:
    # extent 参数严格对应于图像在空间中的实际物理/像素坐标边界
    ax.imshow(background_image, extent=[0, size_x, size_y, 0], aspect='auto')

ax.grid(False)

for group, color in colors.items():
    group_data = spatial_df[spatial_df['leiden'] == group]
    # 进阶技巧：您可以调整 s(点大小) 和 alpha(透明度) 以避免密集区域的过度遮挡
    ax.scatter(group_data['spatial_1'], group_data['spatial_2'], color=color, alpha=0.8, label=group, s=0.5)

# 将图例放置在图像右侧外部，避免遮盖珍贵的组织信息
ax.legend(loc='upper left', bbox_to_anchor=(1.02, 1), markerscale=15, ncol=2, title="Leiden Clusters")
ax.set_title("Custom Spatial Plot overlaid on Tissue Image")
ax.set_xlabel("Spatial X")
ax.set_ylabel("Spatial Y")
plt.gca().set_aspect('equal', adjustable='box')
plt.tight_layout()
plt.show()

8. 添加并展示生物学细胞类型注释

单纯的无监督数字分群（Cluster 0, 1, 2...）在生物学上不够直观。通常我们需要结合已知的 Marker 基因知识，鉴定出这些数字群对应的真实生物学细胞类型（如 T细胞、B细胞、成纤维细胞、肿瘤细胞等）。

假设您已经完成了细胞注释，并将结果保存为 annotation.csv 独立文件（第一列为细胞的 Barcode，并包含名为 Main_CellType 或 Sub_CellType 的注释列）。我们可以轻松地将这份外部注释表合并回数据中，并在 UMAP 和空间组织切片上展示您的鉴定成果。

python

import os

anno_file = "./annotation.csv"
if os.path.exists(anno_file):
    # 读取外部注释文件
    cell_anno = pd.read_csv(anno_file, index_col=0, sep=',')
    
    # 将注释信息无缝合并到 adata 核心对象的 metadata (obs) 中
    cols_to_drop = [col for col in cell_anno.columns if col in adata.obs.columns]
    if cols_to_drop:
        adata.obs.drop(columns=cols_to_drop, inplace=True)
    
    # 将注释信息无缝合并到 adata 核心对象的 metadata (obs) 中
    adata.obs = pd.merge(adata.obs, cell_anno, left_index=True, right_index=True)
    
    # --- 1. 在 UMAP 上展示大类细胞类型 (Main_CellType) ---
    sc.settings.set_figure_params(figsize=(7, 7))
    sc.pl.umap(adata, color=["Main_CellType"], title="UMAP: Main Cell Types")
    sc.pl.umap(adata, color=["Sub_CellType"], title="UMAP: Sub Cell Types")
    
    # --- 2. 在组织空间原位上展示亚群细胞类型 (Sub_CellType) ---
    sc.pl.spatial(adata, color="Main_CellType", library_id=library_id, img_key="DAPI", size=1.5, alpha=1.0, title="Spatial: Main Cell Types")
    sc.pl.spatial(adata, color="Sub_CellType", library_id=library_id, img_key="DAPI", size=1.5, alpha=1.0, title="Spatial: Sub Cell Types")
else:
    print(f"未找到注释文件 '{anno_file}'。\n👉 如果您有注释结果文件，请将其放置在当前目录下并命名为 'annotation.csv'，或者将cell中 annotation.csv 文件名改为您的注释结果的文件名，即可重新运行此单元格查看您的注释图表。")

结语

本教程仅展示了基础的数据读取和原位可视化流程。得益于标准化的 AnnData 对象，您可以无缝接入开源生态内众多强大的空间分析工具（如 Squidpy, Cell2location, Giotto 等），进行空间近邻通讯、细胞互作网络构建、微环境模块发现等高阶深度分析。

单细胞空间转录组： 基础数据分析（基于 Scanpy 的单样本流程） ​

1. 模块简介 ​

2. 输入文件准备 ​

3. 预处理流程：从原始矩阵到 h5ad ​

3.1 读取表达矩阵与空间坐标 ​

3.2 质控 (QC) 与过滤 ​

3.3 归一化、寻找高变基因 (HVG) 与数据标准化 ​

3.4 降维、无监督聚类与差异基因分析 ​

3.5 整合组织切片图像并保存 h5ad ​

4. 读取 h5ad 综合数据对象 (快速开始) ​

5. UMAP 降维与细胞分群展示 ​

6. 原位空间分布展示 ​

7. 叠加组织切片图像展示 ​

7.1 (进阶) 使用 Matplotlib 高度自定义空间图 ​

8. 添加并展示生物学细胞类型注释 ​

单细胞空间转录组：基础数据分析（基于 Scanpy 的单样本流程）

1. 模块简介

2. 输入文件准备

3. 预处理流程：从原始矩阵到 h5ad

3.1 读取表达矩阵与空间坐标

3.2 质控 (QC) 与过滤

3.3 归一化、寻找高变基因 (HVG) 与数据标准化

3.4 降维、无监督聚类与差异基因分析

3.5 整合组织切片图像并保存 h5ad

4. 读取 h5ad 综合数据对象 (快速开始)

5. UMAP 降维与细胞分群展示

6. 原位空间分布展示

7. 叠加组织切片图像展示

7.1 (进阶) 使用 Matplotlib 高度自定义空间图

8. 添加并展示生物学细胞类型注释