从10x Visium到MERFISH:用Scanpy搞定空间转录组数据预处理与可视化的完整流程
从Visium到MERFISH基于Scanpy的空间转录组全流程实战指南空间转录组技术正在彻底改变我们对组织微环境的理解。想象一下你手中握有一张组织切片不仅能看清每个细胞的形态还能精确知道每个位置上哪些基因正在活跃表达——这正是空间转录组带给生物医学研究的革命性视角。不同于传统单细胞测序丢失空间信息的局限这项技术保留了基因表达的空间坐标让研究者能够绘制出组织中分子活动的地理地图。对于刚接触这一领域的生物信息分析师来说最大的挑战往往不是理论理解而是如何快速上手处理这些特殊的数据。本文将带你使用Python生态中最强大的单细胞分析工具Scanpy逐步解析10x Visium和MERFISH两种主流空间转录组数据的处理全流程。无论你是准备分析自己的实验数据还是想要复现文献中的精彩发现这篇实战指南都将成为你实验室工作台上的必备参考手册。1. 实验设计与数据准备空间转录组实验的成功始于明智的技术选择。目前主流平台可分为两大类基于测序的Visium技术和基于成像的MERFISH方法。Visium通过特殊设计的载玻片捕获mRNA每个spot直径55微米通常包含1-10个细胞而MERFISH通过多重荧光原位杂交实现单分子分辨率能精确定位到亚细胞水平。数据获取途径对比技术类型数据来源典型数据量获取方式10x Visium10x Genomics官网5-10GB/样本sc.datasets.visium_sge()MERFISH文献补充材料1-2GB/样本手动下载CSV/Excel对于Visium数据Scanpy提供了便捷的接口直接获取公开数据集。例如获取人类淋巴结数据import scanpy as sc adata sc.datasets.visium_sge(sample_idV1_Human_Lymph_Node) adata.var_names_make_unique()MERFISH数据通常需要从研究论文的补充材料中手动下载。以2019年Xia等人的小鼠下丘脑数据为例import pandas as pd coordinates pd.read_excel(pnas.1912459116.sd15.xlsx, index_col0) counts sc.read_csv(pnas.1912459116.sd12.csv).transpose() adata_merfish counts[coordinates.index, :] adata_merfish.obsm[spatial] coordinates.to_numpy()提示处理MERFISH数据时需特别注意坐标系的统一某些数据集可能使用微米为单位而另一些使用像素坐标。2. 数据质控与预处理实战空间转录组数据的质量控制需要兼顾分子生物学特征和空间特性。与传统单细胞分析不同我们不仅要关注细胞层面的质控指标还需要考虑空间位置上的异常模式。关键质控指标解析线粒体基因比例高于20%可能提示细胞应激或死亡总UMI计数Visium数据通常在5,000-35,000之间检测基因数健康组织一般在1,000-4,000个基因空间离群点检查边缘或孤立区域的高表达点计算质控指标的实用代码# 标记线粒体基因 adata.var[mt] adata.var_names.str.startswith(MT-) sc.pp.calculate_qc_metrics(adata, qc_vars[mt], inplaceTrue) # 可视化质控指标 import seaborn as sns fig, axs plt.subplots(1, 2, figsize(12, 4)) sns.distplot(adata.obs[total_counts], kdeFalse, axaxs[0]) sns.distplot(adata.obs[pct_counts_mt], kdeFalse, axaxs[1])执行过滤的黄金标准# 细胞过滤 sc.pp.filter_cells(adata, min_counts5000) sc.pp.filter_cells(adata, max_counts35000) adata adata[adata.obs[pct_counts_mt] 20] # 基因过滤 sc.pp.filter_genes(adata, min_cells10)空间特异性质控往往被忽视。一个实用的技巧是检查空间坐标上的表达量分布sc.pl.spatial(adata, color[total_counts, n_genes_by_counts], alpha0.7, size1.2)3. 标准化与特征选择策略空间转录组数据的标准化需要特殊考虑。Visium数据由于捕获效率差异spot间的技术变异更为明显而MERFISH数据由于是靶向检测通常更干净但可能受图像分析算法影响。标准化方法对比表方法适用场景优点缺点CPM初步探索计算快未考虑技术变异文库大小校正Visium数据考虑捕获效率忽略基因特性负二项回归复杂实验设计建模精确计算量大Scanpy标准工作流# 基础标准化 sc.pp.normalize_total(adata, target_sum1e4) sc.pp.log1p(adata) # 高变基因选择 sc.pp.highly_variable_genes(adata, flavorseurat, n_top_genes2000) adata adata[:, adata.var.highly_variable]对于空间数据特别推荐检测空间变异基因SVGs。这需要使用SpatialDE等专用工具# 安装空间分析工具 !pip install SpatialDE # 准备输入数据 counts pd.DataFrame(adata.X.todense(), columnsadata.var_names, indexadata.obs_names) coord pd.DataFrame(adata.obsm[spatial], columns[x, y], indexadata.obs_names) # 运行空间差异分析 import SpatialDE results SpatialDE.run(coord, counts) top_svgs results.sort_values(qval).head(50)[g]4. 降维聚类与空间模式解析空间转录组分析的核心价值在于将分子特征映射回组织结构。与传统单细胞分析相比我们需要在降维聚类中特别关注空间连续性。标准分析流程七步法PCA降维sc.pp.pca(adata)邻域图构建sc.pp.neighbors(adata)UMAP可视化sc.tl.umap(adata)Leiden聚类sc.tl.leiden(adata)标记基因鉴定sc.tl.rank_genes_groups()空间可视化sc.pl.spatial()功能注释基于标记基因的GO分析执行基础分析的完整代码# 标准降维聚类 sc.pp.pca(adata) sc.pp.neighbors(adata) sc.tl.umap(adata) sc.tl.leiden(adata, resolution0.6) # 可视化 sc.pl.umap(adata, color[leiden, total_counts])空间特异性的分析技巧是检查聚类结果的空间分布sc.pl.spatial(adata, img_keyhires, colorleiden, size1.5, alpha0.8)对于感兴趣的空间模式可以进一步进行区域特异性分析# 选择特定空间区域 x_min, x_max 3000, 5000 y_min, y_max 2000, 4000 region_mask (adata.obsm[spatial][:,0] x_min) \ (adata.obsm[spatial][:,0] x_max) \ (adata.obsm[spatial][:,1] y_min) \ (adata.obsm[spatial][:,1] y_max) adata_region adata[region_mask].copy() # 区域特异性聚类 sc.tl.leiden(adata_region, key_addedsub_clusters)5. MERFISH数据处理专项技巧MERFISH数据虽然与Visium同属空间转录组但在处理上有其独特之处。由于是成像技术数据通常更稀疏但分辨率更高需要调整分析方法。MERFISH数据处理四步优化归一化策略采用CPM而非文库大小校正特征选择优先使用已知标记基因降维参数减少PCA组分数目10-15聚类分辨率调低Leiden分辨率参数典型MERFISH分析流程# 专用预处理 sc.pp.normalize_per_cell(adata_merfish, counts_per_cell_after1e6) sc.pp.log1p(adata_merfish) # 适度降维 sc.pp.pca(adata_merfish, n_comps15) sc.pp.neighbors(adata_merfish) # 保守聚类 sc.tl.leiden(adata_merfish, resolution0.3) # 空间可视化 sc.pl.embedding(adata_merfish, basisspatial, colorleiden, paletteSet1)MERFISH数据特别适合分析细胞间相互作用。一个实用技巧是计算邻近细胞类型的共现频率from sklearn.neighbors import NearestNeighbors # 计算邻近细胞类型 nbrs NearestNeighbors(n_neighbors5).fit(adata_merfish.obsm[spatial]) distances, indices nbrs.kneighbors(adata_merfish.obsm[spatial]) # 统计邻近关系 cluster_pairs [] for i in range(len(indices)): for j in indices[i][1:]: # 跳过自身 cluster_pairs.append( (adata_merfish.obs.leiden.iloc[i], adata_merfish.obs.leiden.iloc[j]) ) # 生成共现矩阵 co_occurrence pd.crosstab( pd.Series([x[0] for x in cluster_pairs]), pd.Series([x[1] for x in cluster_pairs]) )6. 高级可视化与结果解读优秀的可视化能让空间转录组数据自己讲故事。除了标准绘图外Scanpy支持多种高级展示技巧。空间可视化工具箱函数用途关键参数sc.pl.spatial基础空间图color,size,alphasc.pl.embedding自定义坐标basis,colorsc.pl.heatmap标记基因var_names,groupbysc.pl.dotplot基因模块var_names,groupby展示空间梯度模式的技巧# 创建空间坐标衍生特征 adata.obs[x_coord] adata.obsm[spatial][:,0] adata.obs[y_coord] adata.obsm[spatial][:,1] # 检查基因表达的空间趋势 sc.pl.spatial(adata, color[x_coord, y_coord, Pcp4], ncols3, size1.3)对于多切片分析可以使用scanpy的AnnData拼接功能# 合并多个样本 import anndata adatas [adata_sample1, adata_sample2] adata_combined anndata.concat(adatas, labelsample, keys[sample1, sample2]) # 批次校正 sc.pp.combat(adata_combined, keysample)交互式可视化可以极大提升探索效率。虽然Scanpy本身不支持但可以导出到Napari# 导出为空间数据框架 spatial_df pd.DataFrame(adata.obsm[spatial], columns[x, y], indexadata.obs_names) spatial_df pd.concat([spatial_df, adata.obs], axis1) # 保存为CSV spatial_df.to_csv(spatial_coordinates_with_metadata.csv)7. 疑难解答与性能优化实际分析中常会遇到各种技术挑战。以下是几个常见问题的解决方案问题1Visium数据质控后细胞数过少可能原因过滤阈值过严解决方案检查原始数据质量调整阈值# 动态阈值确定 lower_bound np.percentile(adata.obs[total_counts], 5) upper_bound np.percentile(adata.obs[total_counts], 95)问题2MERFISH聚类结果过于碎片化可能原因分辨率参数过高解决方案尝试不同分辨率for res in [0.1, 0.3, 0.5]: sc.tl.leiden(adata, resolutionres, key_addedfclusters_{res})问题3空间DE分析内存不足可能原因基因数过多解决方案预过滤低表达基因# 保留在至少10%细胞中表达的基因 sc.pp.filter_genes(adata, min_cellsadata.n_obs*0.1)对于大数据集可以采用近似算法提升性能# 使用近似PCA sc.pp.pca(adata, svd_solverarpack) # 降低邻域图大小 sc.pp.neighbors(adata, n_neighbors15)内存优化技巧包括使用稀疏矩阵和分块处理# 转换为稀疏矩阵 from scipy import sparse adata.X sparse.csr_matrix(adata.X) # 分块处理大型数据 chunk_size 1000 for i in range(0, adata.n_obs, chunk_size): chunk adata[i:ichunk_size] process_chunk(chunk)