单细胞空间转录组全流程解析：Visium & MERFISH数据处理全代码实战

本文介绍了单细胞空间转录组数据（Visium和MERFISH）的完整分析流程，包括数据预处理、质量控制、降维聚类、空间可视化和标记基因识别。使用Scanpy、Pandas、Matplotlib和Seaborn等工具，详细展示了从数据加载到结果可视化的全过程。Visium分析涵盖质量控制、标准化、聚类和空间可视化；MERFISH分析则侧重数据加载、归一化和空间聚类。流程提供了完整的Python代码实

2301_76444133

1078人浏览 · 2025-07-24 21:26:47

2301_76444133 · 2025-07-24 21:26:47 发布

单细胞空间转录组数据（Visium和MERFISH）的完整分析流程，涉及数据预处理、质量控制、降维聚类、空间可视化和标记基因识别。

一、环境设置与库导入

import scanpy as sc
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

Scanpy：单细胞数据分析的核心库，支持标准化、降维、聚类等。
Pandas：用于处理Excel和CSV文件。
Matplotlib/Seaborn：可视化工具，用于绘制质量控制图、UMAP图和空间图。

sc.logging.print_versions()  # 打印Scanpy及依赖库版本信息

确保所有工具版本兼容，便于复现实验。

sc.set_figure_params(facecolor="white", figsize=(8, 8))  # 设置默认图表风格
plt.rcParams['figure.figsize'] = (4, 4)  # 设置Matplotlib默认图表尺寸

统一图表风格，提升可视化一致性。

二、Visium数据分析流程

1. 数据加载与初始化

adata = sc.datasets.visium_sge(sample_id="V1_Human_Lymph_Node")

Visium_sge 是10X Genomics提供的空间转录组数据集，sample_id指定样本（此处为人类淋巴结）。
adata 是AnnData对象，包含基因表达矩阵和空间坐标。

adata.var_names_make_unique()  # 确保基因名唯一
adata.var["mt"] = adata.var_names.str.startswith("MT-")  # 标记线粒体基因
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True)

质量控制（QC）：
- 计算每个细胞的total_counts（总读数）、n_genes_by_counts（表达基因数）和pct_counts_mt（线粒体基因占比）。
- 线粒体基因通常用于评估细胞完整性（如线粒体占比过高可能为坏死细胞）。

2. 质量控制与过滤

# 可视化QC指标分布
sns.histplot(adata.obs["total_counts"], kde=False, ax=axs[0])  # 总读数分布
sns.histplot(adata.obs["n_genes_by_counts"], kde=False, bins=60, ax=axs[2])  # 基因数分布

过滤策略：
- sc.pp.filter_cells(adata, min_counts=5000)：保留总读数>5000的细胞。
- sc.pp.filter_cells(adata, max_counts=35000)：去除总读数>35000的细胞（可能为异常细胞）。
- adata = adata[adata.obs["pct_counts_mt"] < 20].copy()：去除线粒体占比>20%的细胞。
- sc.pp.filter_genes(adata, min_cells=10)：保留至少在10个细胞中表达的基因。

3. 数据标准化与特征选择

sc.pp.normalize_total(adata, inplace=True)  # 总读数归一化（使所有细胞总读数相同）
sc.pp.log1p(adata)  # 对数转换，缓解表达值的长尾分布
sc.pp.highly_variable_genes(adata, flavor="seurat", n_top_genes=2000)

标准化：log1p处理后数据更符合正态分布，便于后续分析。
高变基因选择：使用Seurat方法选择表达方差高的2000个基因，这些基因对后续聚类更敏感。

4. 降维与聚类

sc.pp.pca(adata)  # PCA降维
sc.pp.neighbors(adata)  # 构建KNN图（基于PCA结果）
sc.tl.umap(adata)  # UMAP降维
sc.tl.leiden(adata, key_added="clusters", flavor="igraph")  # Leiden聚类（基于UMAP图）

UMAP：非线性降维，用于可视化高维数据。
Leiden聚类：基于图神经网络的聚类算法，可发现更精细的细胞群。

sc.pl.umap(adata, color=["total_counts", "n_genes_by_counts", "clusters"], wspace=0.4)

可视化：用UMAP图展示细胞的分布（total_counts、n_genes_by_counts）和聚类结果。

5. 空间可视化

sc.pl.spatial(adata, img_key="hires", color=["total_counts", "n_genes_by_counts"])
sc.pl.spatial(adata, img_key="hires", color="clusters", size=1.5)

空间图：img_key="hires"表示使用高分辨率图像作为背景，展示细胞的空间分布。
聚类可视化：用不同颜色标注细胞簇，size=1.5控制点的大小。

sc.pl.spatial(adata, img_key="hires", color="clusters", groups=["5", "9"], crop_coord=[7000, 10000, 0, 6000], alpha=0.5)

局部放大：crop_coord指定图像裁剪区域，alpha=0.5调整透明度，突出显示特定簇。

6. 标记基因识别

sc.tl.rank_genes_groups(adata, "clusters", method="t-test") # 根据簇识别标记基因 

sc.pl.rank_genes_groups_heatmap(adata, groups="9", n_genes=10, groupby="clusters")

标记基因：使用t检验找出在特定簇中显著差异表达的基因。
热图：展示簇9的前10个标记基因的表达模式。

sc.pl.spatial(adata, img_key="hires", color=["CR2", "COL1A2", "SYPL1"], alpha=0.7)

关键基因空间表达：用颜色标注特定基因（如CR2、COL1A2）在空间中的表达分布。

三、MERFISH数据分析流程

1. 数据加载与预处理

coordinates = pd.read_excel("./data/pnas.1912459116.sd15.xlsx", index_col=0)
counts = sc.read_csv("./data/pnas.1912459116.sd12.csv").transpose()
adata_merfish = counts[coordinates.index, :].copy()
adata_merfish.obsm["spatial"] = coordinates.to_numpy()

MERFISH：一种高分辨率的单分子成像技术，数据需从Excel和CSV文件加载。
obsm["spatial"]：存储空间坐标（x,y），用于后续空间图可视化。

2. 标准化与降维

sc.pp.normalize_per_cell(adata_merfish, counts_per_cell_after=1e6)  # 每个细胞归一化到1e6个读数
sc.pp.log1p(adata_merfish)  # 对数转换
sc.pp.pca(adata_merfish, n_comps=15)  # PCA降维到15个成分

归一化方法：不同于Visium的全数据归一化，MERFISH使用normalize_per_cell。

3. 聚类与可视化

sc.tl.leiden(adata_merfish, key_added="clusters", resolution=0.5)  # 聚类分辨率参数
sc.pl.umap(adata_merfish, color="clusters")  # UMAP图
sc.pl.embedding(adata_merfish, basis="spatial", color="clusters")  # 空间图

分辨率参数：resolution=0.5控制聚类的粒度（值越大，簇越细）。
空间图：basis="spatial"直接使用原始坐标进行可视化。

完整的代码在下边


import scanpy as sc
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# 打印版本信息
sc.logging.print_versions()

# 设置全局参数
sc.set_figure_params(facecolor="white", figsize=(8, 8))
sc.settings.verbosity = 3
plt.rcParams['figure.figsize'] = (4, 4)  # 设置默认图表尺寸

# ================== Visium数据分析流程 ==================
# 1. 数据加载与初始化
adata = sc.datasets.visium_sge(sample_id="V1_Human_Lymph_Node")
adata.var_names_make_unique()
adata.var["mt"] = adata.var_names.str.startswith("MT-")  # 标记线粒体基因
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True)

# 2. 质量控制与过滤
# 可视化QC指标分布
fig, axs = plt.subplots(1, 4, figsize=(15, 4))
sns.histplot(adata.obs["total_counts"], kde=False, ax=axs[0])
sns.histplot(adata.obs["total_counts"][adata.obs["total_counts"] < 10000],
             kde=False, bins=40, ax=axs[1])
sns.histplot(adata.obs["n_genes_by_counts"], kde=False, bins=60, ax=axs[2])
sns.histplot(adata.obs["n_genes_by_counts"][adata.obs["n_genes_by_counts"] < 4000],
             kde=False, bins=60, ax=axs[3])
plt.show()

# 执行过滤
sc.pp.filter_cells(adata, min_counts=5000)
sc.pp.filter_cells(adata, max_counts=35000)
adata = adata[adata.obs["pct_counts_mt"] < 20].copy()
sc.pp.filter_genes(adata, min_cells=10)

# 3. 数据标准化与特征选择
sc.pp.normalize_total(adata, inplace=True)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, flavor="seurat", n_top_genes=2000)

# 4. 降维与聚类
sc.pp.pca(adata)
sc.pp.neighbors(adata)
sc.tl.umap(adata)
sc.tl.leiden(adata, key_added="clusters", flavor="igraph", directed=False, n_iterations=2)

# 可视化聚类结果
sc.pl.umap(adata, color=["total_counts", "n_genes_by_counts", "clusters"], wspace=0.4)

# 5. 空间可视化
plt.rcParams["figure.figsize"] = (8, 8)
sc.pl.spatial(adata, img_key="hires", color=["total_counts", "n_genes_by_counts"])
sc.pl.spatial(adata, img_key="hires", color="clusters", size=1.5)

# 特定区域放大
sc.pl.spatial(adata, img_key="hires", color="clusters", groups=["5", "9"],
              crop_coord=[7000, 10000, 0, 6000], alpha=0.5, size=1.3)

# 6. 标记基因识别
sc.tl.rank_genes_groups(adata, "clusters", method="t-test")
sc.pl.rank_genes_groups_heatmap(adata, groups="9", n_genes=10, groupby="clusters")

# 关键基因空间表达
sc.pl.spatial(adata, img_key="hires", color=["clusters", "CR2"])
sc.pl.spatial(adata, img_key="hires", color=["COL1A2", "SYPL1"], alpha=0.7)

# ================== MERFISH数据分析流程 ==================
# 1. 数据加载与预处理
coordinates = pd.read_excel("./data/pnas.1912459116.sd15.xlsx", index_col=0)
counts = sc.read_csv("./data/pnas.1912459116.sd12.csv").transpose()
adata_merfish = counts[coordinates.index, :].copy()
adata_merfish.obsm["spatial"] = coordinates.to_numpy()

# 2. 标准化与降维
sc.pp.normalize_per_cell(adata_merfish, counts_per_cell_after=1e6)
sc.pp.log1p(adata_merfish)
sc.pp.pca(adata_merfish, n_comps=15)
sc.pp.neighbors(adata_merfish)
sc.tl.umap(adata_merfish)

# 3. 聚类与可视化
sc.tl.leiden(adata_merfish, key_added="clusters", resolution=0.5,
             n_iterations=2, flavor="igraph", directed=False)

sc.pl.umap(adata_merfish, color="clusters")
sc.pl.embedding(adata_merfish, basis="spatial", color="clusters")