前几天和大家介绍了发表在Cell正刊上的空间转录组分析方法：真漂亮，这是这个月华大空间技术的第二篇Cell了，看看他们的3D时空建模框架——Spateo。这个方法2022年就挂在了预印本上，直到今年11月才正式见刊，也是历程艰辛。最近几期主要和大家分享下，Spateo的使用和部署。
Spateo 是一款专为 3D 空间转录组学（ST）数据分析设计的创新工具，适用于小鼠胚胎（E9.5 和 E11.5 阶段）及其他物种。它突破了传统 3D ST 分析方法的局限，解决了诸如 大规模切片对齐、复杂 3D 形态建模 和 时空动态解析 等难题，助力从单细胞到胚胎全尺度的多维研究。核心功能包括：
3D重建：通过非刚性对齐、多切片优化、网格修正等方法，高效重建全胚尺度3D结构，优于PASTE、SLAT等现有方法。
多尺度空间分析：定义特定空间轴，研究从单细胞到胚胎水平的生物学现象，揭示关键基因（如Slit2、Dbx1）的表达变化。
细胞通信与信号分析：整合细胞间和细胞内交互，解析关键调控网络（如Wnt、Bmp、Fgf8）。
4D动态建模：预测细胞迁移路径，例如心脏细胞的非对称迁移，与Notch信号及调控因子（如Tbx2、Pitx2）相关。
交互式可视化：开发了Spateo-viewer工具，用于3D数据的可视化和操作。
Spateo支持多种ST数据类型，结合Dynamo可定量预测细胞命运转变的时空动力学，并为研究人员提供了开放源码、详细文档及教程（https://github.com/aristoteleo/spateo-release）。Spateo 实现了从单细胞、组织、器官到整个胚胎的多维度研究，开创了 3D ST 数据分析的新篇章。
Spateo的环境部署也比较简单易上手，不懂如何用conda部署安装环境的朋友，可以参考系统性教程02-Linux 安装Anaconda这一期。
环境部署
mamba create -n spateo python=3.9.0 -yconda activate spateopip install spateo-release scanpy stardist imagecodecs cdlib
Spateo 拥有多种功能模块，是空间转录组分析的强大工具。今天重点介绍其 细胞分割模块，这是空间转录组分析中至关重要的一环，对高分辨率的空间转录组数据尤其重要。空间转录组数据本质上是一组具有坐标信息的 mRNA 表达矩阵，但其中 缺失了细胞信息。为了分析细胞级别的表达特征，我们需要通过细胞分割将像素级数据整合为细胞级数据。以百迈客的空间转录组技术为例，他们通过多次对切片成像，提取细胞坐标，最终完成细胞分割并生成细胞级的 mRNA 矩阵，为下游分析提供了高质量的数据。
Spateo 细胞分割模块的特点
支持图像辅助分割：
借助切片图像提取细胞信息，将数据调整到细胞级别，恢复细胞信息。
无图像分割支持：
即使没有图像数据，Spateo 也能基于 RNA 信号进行细胞分割，提供灵活的解决方案。
基于染色结果进行分割
精化对齐：优化染色图像与 RNA 坐标之间的对齐。
基于 Watershed 方法标记细胞核：使用 Watershed 算法识别并标注单个细胞核。
基于深度学习的 StarDist 方法标记细胞核：应用 StarDist 模型完成细胞核的识别和标注。
（可选）整合标签：使用 Watershed 方法的标签补充 StarDist 标签，针对未重叠的部分进行完善。
（可选）扩展标签至细胞质：将细胞核标签扩展到细胞质区域，实现更全面的细胞标注。
import spateo as st # 导入 Spateo 库，用于空间转录组数据分析。import matplotlib.pyplot as plt # 导入 Matplotlib 库，用于可视化。# 配置 Spateo 使用的线程数量为 8，以提升数据处理速度。st.config.n_threads = 8# 配置 Jupyter Notebook 的图表背景为白色，避免默认灰色背景影响图表显示。%config InlineBackend.print_figure_kwargs = {'facecolor' : "w"}# 配置 Jupyter Notebook 图像显示为高分辨率格式，以便在 Retina 屏幕上获得更清晰的图表。%config InlineBackend.figure_format = 'retina'
加载数据
!wget "https://drive.google.com/uc?export=download&id=1nONOaUy7utvtXQ3ZPx7R3TePq2Oo4JFM" -nc -O SS200000135IL-D1.ssDNA.tif!wget "https://drive.google.com/uc?export=download&id=18sM-5LmxOgt-3kq4ljtq_EdWHjihvPUx" -nc -O SS200000135TL_D1_all_bin1.txt.gz将下载的 UMI 计数数据和细胞核染色图像加载到一个 AnnData 对象中。
在进行细胞分割时，我们将使用一个聚合的计数矩阵，其中 AnnData 的行为空间中的 X 坐标，列为 Y 坐标，矩阵中的每个元素表示在每个 X 和 Y 坐标捕获到的 UMI 总数。
adata = st.io.read_bgi_agg('SS200000135TL_D1_all_bin1.txt.gz', 'SS200000135IL-D1.ssDNA.tif',)adata
显示切片
st.pl.imshow(adata, 'stain')染色图像通常已与 RNA 坐标大致对齐，但可能存在小的偏差。大的对齐误差可能导致 UMI 聚合错误，从而影响细胞分割结果。因此，建议在空间转录组技术提供的初始对齐基础上进一步优化对齐。
Spateo 提供两种对齐策略，此处使用较简单的方法，因为测试表明其对当前样本效果良好。但对于其他样本，另一种方法可能表现更佳。
before = adata.layers['stain'].copy()st.cs.refine_alignment(adata, mode='rigid', transform_layers=['stain'])
fig, axes = plt.subplots(ncols=2, figsize=(8, 4), tight_layout=True)axes[0].imshow(before)st.pl.imshow(adata, 'unspliced', ax=axes[0], alpha=0.6, cmap='Reds', vmax=2, use_scale=False, save_show_or_return='return')axes[0].set_title('before alignment')st.pl.imshow(adata, 'stain', ax=axes[1], use_scale=False, save_show_or_return='return')st.pl.imshow(adata, 'unspliced', ax=axes[1], alpha=0.6, cmap='Reds', vmax=2, use_scale=False, save_show_or_return='return')axes[1].set_title('after alignment')
分水岭算法
Spateo包括一个自定义的基于分水岭的方法，用于对染色图像中的细胞核进行分割和标记。总体而言，它结合了全局和局部阈值处理，首先获得细胞核的掩膜（这称为分割），然后使用分水岭算法分配标签（这称为标记）
st.cs.mask_nuclei_from_stain(adata)st.pl.imshow(adata, 'stain_mask')
添加标签
st.cs.find_peaks_from_mask(adata, 'stain', 7)st.cs.watershed(adata, 'stain', 5, out_layer='watershed_labels')fig, ax = st.pl.imshow(adata, 'stain', save_show_or_return='return')st.pl.imshow(adata, 'watershed_labels', labels=True, alpha=0.5, ax=ax)
基于深度学习方案
Spateo采用了多种现有的深度学习方法来进行荧光细胞核分割，如StarDist、Cellpose和DeepCell。我们发现StarDist在这些方法中表现最为一致，因此在这里使用它。需要注意的是，基于深度学习的细胞染色图像分割方法通常将分割和标记结合在一起（即模型输出最终的细胞标签）。
st.cs.stardist(adata, tilesize=-1, equalize=2.0, out_layer='stardist_labels')fig, ax = st.pl.imshow(adata, 'stain', save_show_or_return='return')st.pl.imshow(adata, 'stardist_labels', labels=True, alpha=0.5, ax=ax)
[可选] 增强标签
尽管分水岭方法和StarDist方法在细胞核分割中表现良好，但它们各自也有一些局限性。分水岭方法由于阈值化的性质，往往导致细胞边界粗糙，而StarDist在密集区域有时难以识别细胞核（导致出现“空洞”）。此外，由于我们应用了直方图均衡化（该参数可以通过设置equalize=-1关闭），有时空白区域中的噪声会被放大，错误地被识别为细胞。我们可以通过用分水岭标签增强StarDist标签来缓解这些问题，具体做法是复制那些与任何StarDist标签不重叠的分水岭标签，并删除那些与任何分水岭标签不重叠的StarDist标签。
st.cs.augment_labels(adata, 'watershed_labels', 'stardist_labels', out_layer='augmented_labels')fig, ax = st.pl.imshow(adata, 'stain', save_show_or_return='return')st.pl.imshow(adata, 'augmented_labels', labels=True, alpha=0.5, ax=ax)
[可选] 扩展标签到细胞质
在前面的部分，我们已经识别并标记了染色图像中的单个细胞核。虽然可以直接使用这些标签生成细胞x基因计数矩阵（实际上是核x基因矩阵），但我们建议通过将细胞核标签扩展到细胞质来将核标签转换为细胞标签。
为此，我们将利用ssDNA染色能够弱化染色细胞质的特点，因此可以通过使用宽松的阈值对图像进行阈值处理，从而识别细胞质区域。
st.cs.mask_cells_from_stain(adata, out_layer='stain_cell_mask')st.cs.watershed(adata, 'stain',mask_layer='stain_cell_mask',markers_layer='augmented_labels',out_layer='cell_labels',)
fig, ax = st.pl.imshow(adata, 'stain', save_show_or_return='return')st.pl.imshow(adata, 'cell_labels', labels=True, alpha=0.5, ax=ax)
标记拓展
st.cs.expand_labels(adata, 'augmented_labels', distance=5, max_area=400)fig, ax = st.pl.imshow(adata, 'stain', save_show_or_return='return')st.pl.imshow(adata, 'augmented_labels_expanded', labels=True, alpha=0.5, ax=ax)
获取细胞x基因矩阵
最后，我们将使用分割结果来获得所需的细胞x基因计数矩阵。但首先，我们将扩展所有细胞标签以减轻RNA扩散的影响。扩展的距离将取决于数据中RNA扩散的程度。
st.cs.expand_labels(adata, 'cell_labels', distance=2, out_layer='cell_labels_expanded')fig, ax = st.pl.imshow(adata, 'stain', save_show_or_return='return')st.pl.imshow(adata, 'cell_labels_expanded', labels=True, alpha=0.5, ax=ax)
cell_adata = st.io.read_bgi('SS200000135TL_D1_all_bin1.txt.gz',segmentation_adata=adata,labels_layer='cell_labels_expanded',)cell_adata
RNA分割
如果只有RNA图像可用，而没有任何染色数据，Spateo可以利用RNA信号来识别单个细胞。
import syssys.path.insert(0,"/home/panhailin/software/source/git_hub/spateo-release/")import spateo as stimport matplotlib.pyplot as pltst.config.n_threads = 8%config InlineBackend.print_figure_kwargs = {'facecolor' : "w"}%config InlineBackend.figure_format = 'retina'
加载数据
!wget "https://drive.google.com/uc?export=download&id=18sM-5LmxOgt-3kq4ljtq_EdWHjihvPUx" -nc -O SS200000135TL_D1_all_bin1.txt.gz将下载的UMI计数和细胞核染色图像加载到AnnData对象中。为了进行细胞分割，我们将使用一个聚合的计数矩阵，其中AnnData的obs和var对应空间的X和Y坐标，矩阵的每个元素包含在每个X和Y坐标处捕获的UMI总数。
adata = st.io.read_bgi_agg('SS200000135TL_D1_all_bin1.txt.gz',gene_agg={'nuclear': ['Malat1', 'Neat1']} # Add a layer for nuclear-localized genes)adata
fig, axes = plt.subplots(ncols=3, figsize=(9, 3), tight_layout=True)st.pl.imshow(adata, 'nuclear', ax=axes[0], vmax=2, save_show_or_return='return')st.pl.imshow(adata, 'unspliced', ax=axes[1], vmax=5, save_show_or_return='return')st.pl.imshow(adata, 'X', ax=axes[2], vmax=10)
使用核定位基因识别细胞核
密度分布：正如我们上面所观察到的，图像中存在高RNA密度和低RNA密度的区域。这需要将图像分成几个密度区域，然后分别对每个区域进行分割。否则，算法可能会失去校准，在RNA密集的区域过于敏感，而在RNA稀疏的区域过于严格。Spateo建议首先使用宽松的策略（即将像素分割成多个RNA密度区域），然后手动合并这些区域。
st.cs.segment_densities(adata, 'nuclear', 50, k=3, dk=3, distance_threshold=3, background=False)st.pl.contours(adata, 'nuclear_bins', scale=0.15)
st.pl.imshow(adata, 'nuclear_bins', labels=True)使用VI+BP方法，但也提供了其他方法。我们提供了几种选择（EM+BP、gOTSU、Moran's I、VI+BP）来获取mask，用户可以选择其中之一。
# EM+BPst.cs.score_and_mask_pixels(adata, 'nuclear', k=5, method='EM+BP',)st.pl.imshow(adata, 'nuclear_mask')
# gOTSUst.cs.score_and_mask_pixels(adata, 'nuclear', k=5, method='gauss',)st.pl.imshow(adata, 'nuclear_mask')
# Moran's Ist.cs.moran.run_moran_and_mask_pixels(adata, 'nuclear', k=21)st.pl.imshow(adata, 'nuclear_mask')
# VI+BP (Variational inference + Belief propagation)st.cs.score_and_mask_pixels(adata, 'nuclear', k=5, method='VI+BP',vi_kwargs=dict(downsample=0.1, seed=0, zero_inflated=True))st.pl.imshow(adata, 'nuclear_mask')
这里，我们以EM+BP为例，因此需要重新运行此选项，以替换之前选项生成的nuclear_mask。
# EM+BPst.cs.score_and_mask_pixels(adata, 'nuclear', k=5, method='EM+BP',)st.pl.imshow(adata, 'nuclear_mask')
st.cs.find_peaks_from_mask(adata, 'nuclear', 7)st.cs.watershed(adata, 'nuclear_distances', 1,mask_layer='nuclear_mask',markers_layer='nuclear_markers',out_layer='nuclear_labels')st.pl.imshow(adata, 'nuclear_labels', labels=True)
使用未拼接RNA识别额外的细胞核
密度分箱
st.cs.segment_densities(adata, 'unspliced', 50, k=3, dk=3, distance_threshold=3, background=False)st.pl.contours(adata, 'unspliced_bins', scale=0.15)
st.pl.imshow(adata, 'unspliced_bins', labels=True)分割
然后，像之前对核基因计数所做的那样，使用未拼接RNA识别细胞核。请注意，函数将自动检测RNA密度分箱并相应地调整算法。此外，我们提供了certain_layer参数给st.pp.segmentation.score_and_mask_pixels，以及seed_layer参数给st.pp.segmentation.label_connected_components，用于指示从核定位基因获得的细胞核标签。内部，这些标签被用来进一步帮助识别真实的细胞核。同样，用户可以选择与之前“细胞核定位基因识别”部分中的“分割”部分相同的一个选项来创建前景掩膜。
st.cs.score_and_mask_pixels(adata, 'unspliced', k=5, method='EM+BP',vi_kwargs=dict(downsample=0.1, seed=0),certain_layer='nuclear_labels')st.pl.imshow(adata, 'unspliced_mask')
标记
与之前使用核定位基因时不同，这里我们知道一些初始标签，想要保留这些标签（在某些情况下，还希望将其扩展以填充上面的掩膜）。此外，使用未拼接标签似乎会过度饱和某些RNA密集的区域，以至于很难辨别这些区域中细胞的边界。因此，我们将不使用分水岭方法（它会尝试“填充”整个掩膜），而是使用st.cs.label_connected_components来限制每个标签可以分配的最大区域。请注意，我们提供了seed_label参数，该参数是我们之前使用核定位基因获得的标签。
st.cs.label_connected_components(adata, 'unspliced', seed_layer='nuclear_labels')st.pl.imshow(adata, 'unspliced_labels', labels=True)
筛选细胞质
st.cs.segment_densities(adata, 'X', 50, k=3, distance_threshold=3, dk=5, background=False)st.pl.contours(adata, 'X_bins', scale=0.15)
st.pl.imshow(adata, 'X_bins', labels=True)细胞分割
st.cs.score_and_mask_pixels(adata, 'X', k=7, method='EM+BP',vi_kwargs=dict(downsample=0.1, seed=0),certain_layer='unspliced_labels')st.pl.imshow(adata, 'X_mask')
添加细胞标签
st.cs.label_connected_components(adata, 'X', seed_layer='unspliced_labels')st.pl.imshow(adata, 'X_labels', labels=True)
获取细胞基因矩阵
st.cs.expand_labels(adata, 'X_labels', distance=2, out_layer='X_labels_expanded')st.pl.imshow(adata, 'X_labels_expanded', labels=True)
保存结果
cell_adata = st.io.read_bgi('SS200000135TL_D1_all_bin1.txt.gz',segmentation_adata=adata,labels_layer='X_labels_expanded',)cell_adata
这次生信的大纲内容进行全面的调整，想了解生信的，跟班的，可以看下面👇这个文章
这次可不是只学单细胞，基本上从基础到多组学、空间、机器学习一条龙全打通了
单细胞数据分析需求的可以看👇这个文章
