空间转录组分析实战:从Seurat到Scanpy的Python流程迁移与对比

📅 2026/7/13 15:50:35
空间转录组分析实战:从Seurat到Scanpy的Python流程迁移与对比
1. 空间转录组分析工具生态概览空间转录组技术正在彻底改变我们对组织微环境的理解能力。这项技术不仅能获取基因表达数据还能保留细胞在原始组织中的空间位置信息为研究细胞间相互作用和组织结构提供了全新视角。在数据分析工具领域R语言的Seurat和Python生态的Scanpy是当前最主流的两个选择。Seurat作为空间转录组分析的先驱工具自2018年推出空间分析模块以来已经成为R语言生态中的黄金标准。它提供了一套完整的分析流程从数据读取到可视化都能一站式解决。我刚开始接触空转数据分析时第一个学会的就是Seurat的Load10X_Spatial()函数这个函数能轻松将10x Genomics的Space Ranger输出转换为可分析的Seurat对象。而Scanpy作为Python生态中的对应方案近年来发展迅猛。它基于AnnData数据结构与Python强大的机器学习生态无缝集成。我在去年一个需要整合深度学习模型的项目中就深有体会——当分析流程需要从Seurat切换到PyTorch时数据转换的麻烦程度让我下定决心掌握Scanpy。这两个工具在核心功能上高度相似都支持数据读取与质量控制标准化与高变基因筛选降维聚类分析空间可视化差异表达分析但它们的实现方式和扩展性各有特色。Seurat在交互式探索和可视化方面更胜一筹而Scanpy在处理大规模数据和与AI框架集成时表现更好。下面这个简单对比表格概括了它们的主要差异特性Seurat (R)Scanpy (Python)数据结构Seurat对象AnnData对象标准化方法SCTransformpp.normalize_total高变基因筛选FindVariableFeaturespp.highly_variable_genes降维RunPCA/RunUMAPpp.pca/tl.umap聚类算法FindNeighbors/FindClusterspp.neighbors/tl.leiden可视化SpatialFeaturePlotsc.pl.spatial深度学习兼容性有限优秀在实际项目中选择哪种工具往往取决于团队的技术栈和分析需求。如果你已经熟悉R语言且不需要复杂的机器学习扩展Seurat可能是更稳妥的选择。但如果你计划进行大规模数据分析或需要与TensorFlow/PyTorch集成Scanpy的优势就会凸显出来。2. 数据格式转换与Scanpy对象创建从Seurat转向Scanpy的第一步就是数据格式的转换。10x Genomics的空间转录组数据通常包含三个关键部分表达矩阵filtered_feature_bc_matrix.h5、空间坐标信息spatial/tissue_positions_list.csv和组织切片图像spatial/tissue_lowres_image.png。在Seurat中我们使用Load10X_Spatial()一键读取而在Scanpy中则需要分步处理。让我用一个真实案例来说明这个过程。最近在分析一个脑胶质瘤样本GSM5833536时我首先在R中用Seurat完成了初步分析library(Seurat) GBM4 - Load10X_Spatial( data.dir GBM4_spaceranger_out, filename filtered_feature_bc_matrix.h5, slice GBM4 )要将这个数据迁移到Python环境我们需要先将Seurat对象转换为Scanpy可读的格式。这里推荐使用rpy2库在Python中调用R函数完成转换import rpy2.robjects as ro from rpy2.robjects import pandas2ri from rpy2.robjects.conversion import localconverter # 激活R-Python转换环境 pandas2ri.activate() # 读取Seurat对象 ro.r(load(GBM4.rdata)) # 假设已保存为RData文件 with localconverter(ro.default_converter pandas2ri.converter): adata ro.r( library(Seurat) library(SeuratDisk) as.SingleCellExperiment(GBM4) )不过更直接的方法是直接从原始数据构建Scanpy对象。以下是标准的Python实现import scanpy as sc import pandas as pd # 读取表达矩阵 adata sc.read_10x_h5(filtered_feature_bc_matrix.h5) # 读取空间坐标 positions pd.read_csv(spatial/tissue_positions_list.csv, headerNone, index_col0) positions.columns [in_tissue, array_row, array_col, pxl_row_in_fullres, pxl_col_in_fullres] # 将坐标信息添加到AnnData对象 adata.obsm[spatial] positions[[pxl_row_in_fullres, pxl_col_in_fullres]].values # 添加组织切片信息可选 import imageio img imageio.imread(spatial/tissue_lowres_image.png) adata.uns[spatial] {GBM4: {images: {lowres: img}}}在这个过程中最容易踩的坑是坐标系统的匹配问题。我发现有些实验室提供的坐标文件可能使用不同的参考系这时需要手动调整坐标方向。一个实用的检查方法是绘制spot分布图import matplotlib.pyplot as plt plt.scatter(adata.obsm[spatial][:,0], -adata.obsm[spatial][:,1], s1) plt.gca().set_aspect(equal) plt.show()3. 标准化与高变基因筛选的对比实现数据标准化是空间转录组分析中最关键的预处理步骤之一。Seurat的SCTransform和Scanpy的标准流程在理念上相似但实现细节有显著差异。在Seurat中我们通常使用SCTransform进行标准化GBM4 - SCTransform(GBM4, assay Spatial, verbose FALSE)这个函数一次性完成了三件事1) 归一化处理2) 方差稳定转换3) 高变基因筛选。而在Scanpy中这三个步骤需要分开执行# 基本标准化 sc.pp.normalize_total(adata, target_sum1e4) sc.pp.log1p(adata) # 高变基因筛选 sc.pp.highly_variable_genes( adata, flavorseurat, n_top_genes2000, subsetTrue )我在实践中发现几个值得注意的差异点基因选择策略Seurat的FindVariableFeatures默认使用vst方法而Scanpy支持更多选择seurat、cell_ranger等。在分析小鼠脑数据时我发现当细胞数超过1万时使用flavorcell_ranger效果更好。标准化强度Scanpy的normalize_total默认将所有细胞的计数归一化到中位数而Seurat的SCTransform会考虑基因-细胞间的复杂关系。对于质量差异大的样本SCTransform可能更稳健。内存消耗处理大型数据集时Scanpy的内存效率通常更高。我曾比较过同一个样本约5000个spot的处理Seurat占用约8GB内存而Scanpy仅需3GB。下面是一个更接近SCTransform效果的Scanpy实现方案# 更接近SCTransform的标准化流程 import numpy as np from scipy.sparse import csr_matrix # 负二项回归归一化 counts adata.X.copy() if isinstance(adata.X, csr_matrix) else csr_matrix(adata.X) size_factors np.array(counts.sum(axis1)).ravel() / np.median(counts.sum(axis1)) norm_counts counts / size_factors[:, None] # 方差稳定转换 (近似值) norm_counts.data np.log1p(norm_counts.data * 1e4) adata.X norm_counts # 基于分散度的高变基因筛选 sc.pp.highly_variable_genes(adata, flavorseurat_v3)在实际项目中我通常会保存两种标准化结果进行比较。特别是当后续分析涉及差异表达时不同标准化方法可能导致显著不同的结果。4. 降维与聚类的Python实现降维和聚类是揭示空间转录组数据生物学模式的核心步骤。Seurat和Scanpy在这一环节的语法非常相似但参数设置和默认算法有些微差别。4.1 PCA降维在Seurat中PCA降维直接作用于SCT assayGBM4 - RunPCA(GBM4, assay SCT, verbose FALSE)Scanpy中的对应操作更加模块化# 首先计算PCA sc.pp.pca(adata, n_comps50, svd_solverarpack) # 可视化方差解释比例 plt.plot(np.cumsum(adata.uns[pca][variance_ratio])) plt.xlabel(PCs) plt.ylabel(Cumulative Variance Explained) plt.axvline(x15, colorred, linestyle--)一个实用技巧是通过肘部法则确定主成分数。我通常绘制累积方差图选择解释大部分变异的最小PC数。对于大多数空间转录组数据15-30个PC已经足够。4.2 邻域图构建与聚类Seurat使用FindNeighbors和FindClusters进行聚类GBM4 - FindNeighbors(GBM4, reduction pca, dims 1:10) GBM4 - FindClusters(GBM4, resolution 0.4)Scanpy的等价实现为# 构建KNN图 sc.pp.neighbors(adata, n_neighbors15, n_pcs10) # Leiden聚类 sc.tl.leiden(adata, resolution0.4, key_addedclusters)这里有几个参数需要特别注意n_neighbors控制局部邻域大小值太大会丢失细节太小会引入噪声。我通常在10-30之间调整。resolution决定聚类粒度值越大簇越多。对于空间数据0.4-1.0范围通常合适。随机种子为确保结果可重复务必设置random_state参数。4.3 UMAP/tSNE降维可视化Seurat中的UMAP实现GBM4 - RunUMAP(GBM4, reduction pca, dims 1:10)Scanpy中的对应方法# UMAP降维 sc.tl.umap(adata, min_dist0.3, spread1.0) # 可视化 sc.pl.umap(adata, color[clusters, n_counts], wspace0.4)在空间转录组分析中我强烈建议同时检查基因表达空间和物理空间的聚类结果。这能帮助识别是真实的生物模式还是技术假象# 空间坐标可视化 sc.pl.spatial( adata, img_keyhires, colorclusters, size1.5, alpha0.7 )5. 空间可视化与标记基因分析空间可视化是空间转录组分析的最大特色。Seurat的SpatialFeaturePlot和Scanpy的sc.pl.spatial都能将基因表达映射回组织切片但它们的定制化程度有所不同。5.1 基础空间可视化在Seurat中绘制标记基因的空间分布SpatialFeaturePlot(GBM4, features c(SOX10, SOD2))Scanpy的等价实现sc.pl.spatial( adata, img_keyhires, color[SOX10, SOD2], alpha0.7, bwTrue # 灰度显示背景 )我发现Scanpy在图像处理上更灵活。例如可以轻松调整图像透明度来突出表达信号# 自定义可视化 with plt.rc_context({figure.figsize: (8, 4)}): fig, axs plt.subplots(1, 2) for i, gene in enumerate([SOX10, SOD2]): sc.pl.spatial( adata, img_keyhires, colorgene, axaxs[i], showFalse, alpha_img0.3, alpha0.9, size1.8, titlegene, bwTrue ) plt.tight_layout()5.2 标记基因识别Seurat使用FindAllMarkers进行差异表达分析markers - FindAllMarkers(GBM4, only.pos TRUE)Scanpy提供了多种差异表达分析方法# 使用t-test sc.tl.rank_genes_groups( adata, clusters, methodt-test, key_addedt_test ) # 使用Wilcoxon检验更接近Seurat默认方法 sc.tl.rank_genes_groups( adata, clusters, methodwilcoxon, key_addedwilcoxon ) # 可视化top标记基因 sc.pl.rank_genes_groups_heatmap( adata, groups5, n_genes10, groupbyclusters, keywilcoxon )对于大型数据集我推荐使用更高效的logreg方法sc.tl.rank_genes_groups( adata, clusters, methodlogreg, n_genes200, key_addedlogreg )5.3 空间共表达模式分析Scanpy的一个独特优势是可以方便地整合空间信息进行共表达分析。例如我们可以计算基因对的空间相关性from scipy.stats import spearmanr def spatial_correlation(adata, gene1, gene2): 计算两个基因的空间表达相关性 idx1 adata.var_names.get_loc(gene1) idx2 adata.var_names.get_loc(gene2) exp1 adata.X[:, idx1].toarray().flatten() exp2 adata.X[:, idx2].toarray().flatten() return spearmanr(exp1, exp2).correlation # 示例计算SOX10和OLIG2的空间相关性 corr spatial_correlation(adata, SOX10, OLIG2) print(f空间相关性{corr:.3f})6. 计算效率与扩展性对比在实际应用中工具的计算效率和扩展性往往是关键考量因素。我从三个方面对比了Seurat和Scanpy的性能表现。6.1 内存占用在处理大型空间数据集时如Visium HD包含数万个spot内存管理尤为关键。在我的测试中对于约5,000个spot的数据集Seurat占用约6-8GB内存Scanpy占用约3-5GB内存对于约20,000个spot的数据集Seurat需要30GB内存Scanpy保持在15-20GB范围这种差异主要源于Python稀疏矩阵处理的优化。Scanpy默认使用scipy.sparse格式存储数据而Seurat在R4.0之后才完全支持稀疏矩阵。6.2 计算速度下表比较了关键步骤在相同硬件16核CPU32GB内存上的运行时间步骤Seurat (s)Scanpy (s)数据读取4528标准化12065PCA降维8542聚类 (5,000 cells)3018UMAP可视化2512Scanpy在大多数步骤上都有2倍左右的优势特别是在PCA和UMAP计算上。这是因为其底层调用了更优化的数值计算库。6.3 与深度学习框架集成这是Scanpy最显著的优势。例如我们可以直接将预处理好的数据输入PyTorch模型import torch from torch.utils.data import Dataset class SpatialDataset(Dataset): def __init__(self, adata): self.X torch.FloatTensor(adata.X.toarray()) self.positions torch.FloatTensor(adata.obsm[spatial]) def __len__(self): return len(self.X) def __getitem__(self, idx): return self.X[idx], self.positions[idx] dataset SpatialDataset(adata) loader torch.utils.data.DataLoader(dataset, batch_size32)这种无缝集成使得复杂分析流程如图神经网络应用于空间数据成为可能。相比之下Seurat需要通过Reticulate等桥接工具才能调用Python函数增加了复杂度。7. 实战经验与常见问题解决在长期使用这两种工具的过程中我积累了一些宝贵的实战经验特别是针对迁移过程中的常见问题。7.1 数据转换的陷阱当需要将Seurat对象转换为Scanpy对象时直接转换基因名可能会出问题。因为R的Make.unique和Python的索引处理不同基因名中的特殊字符如-在两者中的处理方式不同一个稳健的解决方案是# 从Seurat导出为h5ad ro.r( library(SeuratDisk) SaveH5Seurat(GBM4, filename GBM4.h5Seurat) Convert(GBM4.h5Seurat, dest h5ad) ) # 在Python中读取 adata sc.read_h5ad(GBM4.h5ad)7.2 批次效应处理当整合多个样本时Scanpy的批次校正方法更丰富# 使用BBKNN进行批次校正 sc.external.pp.bbknn(adata, batch_keysample) # 或者使用Harmony import harmonypy adata.obsm[X_pca_harmony] harmonypy.run_harmony( adata.obsm[X_pca], adata.obs, sample )相比之下Seurat的IntegrateData在空间数据上有时会过度校正丢失真实的空间变异信号。7.3 空间特异性分析Scanpy可以方便地结合空间信息进行特异性分析。例如计算空间自相关from esda.moran import Moran from libpysal.weights import DistanceBand # 创建空间权重矩阵 coords adata.obsm[spatial] w DistanceBand(coords, threshold300) # 计算SOX10的空间自相关 moran Moran(adata[:, SOX10].X.toarray().flatten(), w) print(fMorans I: {moran.I:.3f}, p-value: {moran.p_norm:.4e})这种分析在Seurat中需要额外安装R包并编写复杂代码。7.4 可视化定制Scanpy的绘图基于matplotlib可以深度定制。例如创建复合空间图fig, ax plt.subplots(1, 2, figsize(12, 6)) # 左侧聚类结果 sc.pl.spatial( adata, img_keyhires, colorclusters, axax[0], showFalse, titleClusters ) # 右侧标记基因表达 sc.pl.spatial( adata, img_keyhires, colorSOX10, axax[1], showFalse, titleSOX10 Expression, color_mapmagma ) plt.tight_layout() plt.savefig(combined_spatial.pdf, dpi300)这种灵活的可视化组合在发表级图片制作中非常有用。