空间转录组(ST,Spatial transcriptomics)技术正在革新探索组织空间结构的方式。目前,ST数据分析通常仅限于单个2D组织切片,限制了我们理解3D空间中发生的生物过程。在这里,作者提出了STitch3D,一个统一的框架,联合多个ST切片重建3D生物结构。通过对多个切片进行联合建模,将其与单细胞RNA测序数据进行整合,STitch3D同时识别出具有一致基因表达水平的3D空间区域,并揭示3D细胞类型分布。STitch3D能从批次效应中区分切片之间的生物异质性,并有效地借用切片之间的信息来组成3D模型。实验展示了STitch3D在构建综合3D对象方面的性能,可以在整个组织区域甚至整个生物体中进行3D分析。STitch3D的输出可用于多个下游任务,从而全面了解生物系统。
来自:Construction of a 3D whole organism spatial atlas by joint modelling of multiple slices with deep neural networks
工程地址:https://github.com/YangLabHKUST/STitch3D
目录
- 背景概述
- 方法
- 方法概述
- 预处理
- 学习空间信息
- 基因的选择及细胞特异性基因表达谱矩阵的构建
- 网络结构
- 模型训练
- 基因表达的去噪或插补
- 评估
- 在真实数据上评估
- 在模拟数据上评估
- 基准研究
- 实验
- 重建成年小鼠大脑
背景概述
ST技术可实现完整组织中高通量基因表达谱分析。最近,包含组织中多个平行2D切片的数据集已经获得积累,但大多数现有工具仅提供单个2D切片的分析。器官在三维空间中被组织成复杂的结构,而生物过程很少发生在单一的二维平面上。例如,在胚胎发育过程中,形态的变化在3D空间中扩散。仅2D角度极大地限制了我们对生物信号和过程的解释。
切片描绘了一个2D的转录组景观,多个切片的联合建模为描绘生物系统的三维场景提供了机会。为此,需要对多个切片进行联合分析,以重建详细的3D组织结构。具体来说,以下两个3D分析任务是基本的。首先是识别具有相似基因表达spot的3D空间区域,以揭示组织结构。结果可以进行下游分析,比如通过3D的patterns检测区域相关的基因。第二个任务是通过整合多个ST切片和scRNA-seq图谱来推断3D细粒度细胞类型分布。
现有的ST技术在每个spot上通常包含多个细胞,导致分辨率相对较低。利用来自scRNA-seq图谱的信息,3D cell-type deconvolution任务在空间点上分解细胞混合物,从而实现更高的3D重建分辨率,从而深入了解特定细胞型富集区域的生物功能。
在这篇文章中,作者提出了使用多个切片来表征3D组织结构的STitch3D。它在一个统一的框架中解决了上述两个3D分析任务。通过对多个切片的基因表达和空间位置进行有效建模,STitch3D可以从批次效应中区分切片之间的生物异质性,并整合跨切片的信息以组建逼真的3D组织。基准研究和真实数据示例说明了STitch3D最先进的准确性和跨多个切片借用信息的能力。现有的ST工具大多是通过执行spatial domain检测任务或cell-type deconvolution任务来分析单个切片,与之相比,STitch3D可以通过精确的3D domain检测和3D cell-type deconvolution结果来构建细胞组织结构的3D视图。使用不同的数据集,都展示了STitch3D在组织规模甚至整个生物体规模的3D重建能力。
- 空间域检测:3D空间下的spot细胞类型识别或者聚类(基于表征的判别式)
- 细胞型反卷积:去噪或者插补(基于表征的生成式)
方法
方法概述
STitch3D是一种基于深度学习的方法,它使用多个2D ST切片来重建3D组织结构(图1)。STitch3D的输入是多个ST切片和匹配的scRNA-seq。预处理步骤包括通过对切片建立三维空间点坐标和构建全局三维邻域图。在这些步骤之后,STitch3D被训练来整合所有切片的信息。引入共享潜在空间来提取有意义的生物异质性并促进批次效应的去除。在embedding空间中,每个点都有自己的表示,用来共同完成空间域识别和细胞型反卷积任务。
具体来说,STitch3D通过利用spot的3D邻域图,将多个切片中spot的基因表达和空间信息映射到共享的潜在空间中。引入另一个神经网络从潜在表征中推断细胞类型。STitch3D可以通过结合估计的细胞类型和细胞类型特异性基因表达谱来重建ST基因表达。训练后,STitch3D输出空间点表示和细胞类型。这些表示用于聚类算法的三维空间域识别,而细胞类型群体有助于恢复三维细胞类型分布。有了这些输出,STitch3D可以进行各种下游分析。
- 图1a:来自多个ST组织切片的原始数据,以及来自参考scRNA-seq数据集的细胞类型特异性基因表达谱,作为STitch3D的输入。
- 图1b:STitch3D的预处理步骤包括:对来自不同组织切片的spot进行对齐,构建spot的三维位置,构建全局三维graph。STitch3D的主模型结合这些结构进行三维空间域识别和三维细胞型反卷积的表示学习。
- 图1c:STitch3D输出三维空间区域识别结果和组织中不同细胞类型的三维空间分布估计。STitch3D还支持多种下游分析,包括空间轨迹推断、低质量基因表达测量的去噪、虚拟组织切片的生成以及具有3D空间表达模式的基因识别。
- 图1d:STitch3D联合对多个切片进行建模,并利用基于图注意力的神经网络,利用三维空间信息学习spot和细胞类型的潜在表征。
预处理
STitch3D首先结合了迭代最近点算法(ICP)或PASTE算法,通过空间点的配准来对齐多个切片。STitch3D通过空间点的注册来对齐多个切片。默认情况下,STitch3D采用ICP算法对切片进行两两对齐。在对齐之前,STitch3D根据邻域的数量检测每个切片中的边缘点。例如,对于使用Visium平台获取的数据集,一个非边缘点应该有六个邻近点,为了保证鲁棒性,定义了一个边缘点,其邻域点的个数小于5且大于1。在用ST平台分析的数据集中,一个非边缘点应该有8个邻居点。相应地,如果spot具有少于7个且多于1个邻点,则归类为边缘点。接下来,假设边缘点通常位于组织的边界上,STitch3D使用ICP算法对两个切片中的边缘点进行对齐。令 P = { p 1 , . . , p m } P=\left\{p_{1},..,p_{m}\right\} P={p1,..,pm}为源切片中边缘点位置的集合, Q = { q 1 , . . . , q k } Q=\left\{q_{1},...,q_{k}\right\} Q={q1,...,qk}为目标切片中边缘点位置的集合,其中, p i , q i ∈ R 2 p_{i},q_{i}\in \mathbb{R}^{2} pi,qi∈R2。ICP算法通过迭代执行以下步骤来对齐点,直到收敛:
- 对于 P P P中的每个点,在目标点云 Q Q Q中找到距离其最近的点,形成集合 Q ′ Q' Q′: { q 1 ′ , . . . , q m ′ } \left\{q_{1}',...,q_{m}'\right\} {q1′,...,qm′}
- 通过求解最小二乘问题求出旋转和平移最优的变换: R ^ , v ^ = a r g m i n R ∈ R 2 × 2 , v ∈ R 2 ∑ i = 1 m ∣ ∣ ( R p i + v ) − q i ′ ∣ ∣ \widehat{R},\widehat{v}=argmin_{R\in\mathbb{R}^{2\times 2},v\in\mathbb{R}^{2}}\sum_{i=1}^{m}||(Rp_{i}+v)-q'_{i}|| R ,v =argminR∈R2×2,v∈R2i=1∑m∣∣(Rpi+v)−qi′∣∣其中约束 R ⊤ R = I R^{\top}R=I R⊤R=I
- 应用 R ^ p i + v ^ \widehat{R}p_{i}+\widehat{v} R pi+v 到源点云 P P P上
对于多个切片,STitch3D按顺序使用ICP算法进行两两对齐。在多个切片对齐后,STitch3D沿着新引入的z轴(描述切片之间的距离)组装2D切片,以构建点的三维空间位置。利用点的三维空间位置,STitch3D为所有切片中的点构建一个全局三维邻域图。默认情况下,如果两个点的距离小于一个切片中最近的两个点之间距离的1.1倍,则在全局3D graph中生成一条边相连。
然后,使用STitch3D模型进行3D空间区域检测和3D细胞分型反卷积。STitch3D通过整合ST数据集的多个切片和scRNA-seq参考数据集的信息来执行3D分析。在ST数据中,使用 s = 1 , 2 , . . . , S s=1,2,...,S s=1,2,...,S作为组织切片的索引。对于切片 s s s,记spot数为 N s N_{s} Ns,观测到的基因表达矩阵为 T s = [ Y n , g s ] ∈ R N s × G T^{s}=[Y^{s}_{n,g}]\in\mathbb{R}^{N_{s}\times G} Ts=[Yn,gs]∈RNs×G, n = 1 , 2 , . . . , N s n=1,2,...,N_{s} n=1,2,...,Ns为细胞, g = 1 , 2 , . . . , G g=1,2,...,G g=1,2,...,G为基因。在注释的scRNA-seq reference中,将获得的细胞类型特定的基因表达矩阵记为 V = [ V c , g ] ∈ R C × G V=[V_{c,g}]\in\mathbb{R}^{C\times G} V=[Vc,g]∈RC×G, c = 1 , 2 , . . . , C c=1,2,...,C c=1,2,...,C代表细胞类型的index。其中的向量 v c ∈ R G v_{c}\in\mathbb{R}^{G} vc∈RG表示类型 c c c细胞的平均基因表达。
学习空间信息
为了整合所有组织切片的基因表达水平和空间位置进行3D分析,STitch3D首先将spot的基因表达信息编码到共享的潜在空间中。为了清晰,我们将各个切片上的点的index s = 1 , 2 , . . . , S s=1,2,...,S s=1,2,...,S连接起来,构建全局index i = 1 , 2 , . . . , N i=1,2,...,N i=1,2,...,N,其中 N = N 1 + ⋅ ⋅ ⋅ + N S N=N_{1}+\cdot\cdot\cdot+N_{S} N=N1+⋅⋅⋅+NS。记 A = [ A i , j ] ∈ R N × N A=[A_{i,j}]\in\mathbb{R}^{N\times N} A=[Ai,j]∈RN×N为3D全局邻居graph, A i , j A_{i,j} Ai,j取值为0或1。令 Y = [ Y i , g ] ∈ R N × G Y=[Y_{i,g}]\in\mathbb{R}^{N\times G} Y=[Yi,g]∈RN×G为所有切片上连接的基因表达计数矩阵。具体地说,spot的潜在表示为: X i , g = l o g ( Y i , g ∑ g = 1 G Y i , g × 1 0 4 + 1 ) Z = f Z ( X , A ) X_{i,g}=log(\frac{Y_{i,g}}{\sum_{g=1}^{G}Y_{i,g}}\times10^{4}+1)\\Z=f_{Z}(X,A) Xi,g=log(∑g=1GYi,gYi,g×104+1)Z=fZ(X,A)首先对计数矩阵 Y Y Y进行归一化和对数变换来保证网络的稳定性, X ∈ R N × G X\in\mathbb{R}^{N\times G} X∈RN×G被编码到联合空间下的表示 Z ∈ R N × p Z\in\mathbb{R}^{N\times p} Z∈RN×p。然后,STitch3D根据spot的潜在表示生成细胞类型比例,记spot i i i中 c c c类型细胞的比例为 β i , c \beta_{i,c} βi,c。记 β i = [ β i , 1 , β i , 2 , . . . , β i , C ] ⊤ ∈ R C \beta_{i}=[\beta_{i,1},\beta_{i,2},...,\beta_{i,C}]^{\top}\in\mathbb{R}^{C} βi=[βi,1,βi,2,...,βi,C]⊤∈RC表示spot i i i的细胞类型比例, Z i ∈ R p Z_{i}\in\mathbb{R}^{p} Zi∈Rp为spot i i i的embedding, i = 1 , 2 , . . . , N i=1,2,...,N i=1,2,...,N。我们假设细胞类型比例 β i = f β ( Z i ) \beta_{i}=f_{\beta}(Z_{i}) βi=fβ(Zi)。值得注意的是,通过GAT f Z ( ⋅ ) f_{Z}(\cdot) fZ(⋅)得到的 Z Z Z是包含3D信息的。因此, β i \beta_{i} βi也纳入了空间信息。
为了考虑ST数据与scRNA-seq数据之间的技术噪声批次效应,以及多个ST切片之间的批次效应,STitch3D进一步引入 α i s ∈ R \alpha_{i}^{s}\in\mathbb{R} αis∈R和 γ g s ∈ R \gamma_{g}^{s}\in\mathbb{R} γgs∈R,前者表示slice-spot特异性效应,后者表示slice-gene特异性效应。结合细胞类型特异性基因表达谱矩阵 V = [ V c , g ] V = [V_{c, g}] V=[Vc,g],细胞类型比例 β i = [ β i , c ] β_{i} = [β_{i, c}] βi=[βi,c],以及 α i s α^s_i αis和 γ g s γ^s_g γgs两种效应,STitch3D可以用以下模型重建ST数据中观察到的计数: Y i , g ∼ P o i s s o n ( l i , λ i , g ) λ i , g = e x p [ l o g ( ∑ c = 1 C β i , c V c , g ) + α i s + γ g s ] Y_{i,g}\sim Poisson(l_{i},\lambda_{i,g})\\ \lambda_{i,g}=exp[log(\sum_{c=1}^{C}\beta_{i,c}V_{c,g})+α^s_i+γ^s_g] Yi,g∼Poisson(li,λi,g)λi,g=exp[log(c=1∑Cβi,cVc,g)+αis+γgs]其中 l i l_i li是在spot i i i观察到的转录数。我们也可以采用负二项分布NB来模拟 γ i , g \gamma_{i,g} γi,g作为替代。STitch3D将泊松分布设置为默认。根据经验,作者发现基于泊松分布和NB分布的STitch3D模型具有相当的性能。
作者利用slice特异性网络从共享空间生成slice-spot特异性效应 α i s α^s_i αis: α i s = f α ( Z i , s ) \alpha_{i}^{s}=f_{\alpha}(Z_{i},s) αis=fα(Zi,s)其中 s s s为slice的label,我们建模 α i s \alpha_{i}^{s} αis如下:
- 假设slice-spot特异性效应 α i s \alpha_{i}^{s} αis与基因表达信息(比如 Z i Z_i Zi)和slice label信息(在 s s s中)都有关;
- slice-gene特异性效应 γ g s \gamma_{g}^{s} γgs建模为可训练参数;
损失函数为: L = − 1 N ∑ i = 1 N ∑ g = 1 G [ Y i , g l o g ( l i , λ i , g ) − l i λ i , g ] L=-\frac{1}{N}\sum_{i=1}^{N}\sum_{g=1}^{G}[Y_{i,g}log(l_{i},\lambda_{i,g})-l_{i}\lambda_{i,g}] L=−N1i=1∑Ng=1∑G[Yi,glog(li,λi,g)−liλi,g]在STitch3D中,作者考虑了一个正则化器来鼓励保存多个切片上的生物异质性。具体来说,正则化器设计为: R A E = 1 N ∑ i = 1 N ∣ ∣ f X ( Z i , s ) − X ∣ ∣ 2 R_{AE}=\frac{1}{N}\sum_{i=1}^{N}||f_{X}(Z_{i},s)-X||_{2} RAE=N1i=1∑N∣∣fX(Zi,s)−X∣∣2其中,网络 f X ( ⋅ , s ) f_{X}(\cdot,s) fX(⋅,s)用于重建 X i X_{i} Xi。两个网络 f Z ( ⋅ ) f_{Z}(\cdot) fZ(⋅)和 f X ( ⋅ , s ) f_{X}(\cdot,s) fX(⋅,s)共同构成 X X X和 Z Z Z之间的自编码器。在正则化器 R A E R_{AE} RAE中,网络 f X ( ⋅ , s ) f_{X}(⋅,s) fX(⋅,s)被设计用于考虑切片之间的批次效应,从而鼓励网络 f Z ( ⋅ ) f_{Z}(⋅) fZ(⋅)提取生物信息。STitch3D中的最终目标函数由 L + k A E R A E L + k_{AE}R_{AE} L+kAERAE给出,其中 k A E k_{AE} kAE为正则化器 R A E R_{AE} RAE的系数,固定为0.1。
通过上述设计,STitch3D应用SGD来更新模型中的可训练参数。训练后,STitch3D输出表征 Z i Z_i Zi和细胞类型比例 β i β_i βi。利用如Louvain算法和基于高斯混合模型的聚类算法,完成三维空间域识别任务。结合获得的三维空间位置,这些识别的簇揭示了所有切片的三维空间区域。同时,细胞类型比例 β i β_i βi显示了不同细胞类型的三维空间分布,为更高空间分辨率的ST研究提供了一个全面的视角。
基因的选择及细胞特异性基因表达谱矩阵的构建
STitch3D根据scRNA-seq参考数据选择基因。具体来说,对于注释的scRNA-seq中的每种细胞类型数据,STitch3D使用t-test来查找默认设置 K = 500 K = 500 K=500的top-K标记基因。然后将所有细胞类型的top标记基因连接起来,以获得用于STitch3D表示信息基因的完整列表。
STitch3D从 C C C个细胞类型中选择 G G G个基因,构建细胞类型特异性基因表达谱矩阵 V ∈ R C × G V\in\mathbb{R}^{C×G} V∈RC×G。首先,STitch3D对参考scRNA-seq数据集中每个细胞的基因表达量进行归一化处理,使每个细胞的基因表达量之和等于1。然后,计算 c c c细胞类型的行向量 V c V_c Vc,即 c c c细胞类型的平均表达特征。如果参考scRNA-seq数据集中包含不同批次的细胞,则STitch3D为每批次计算一个细胞类型特异性基因表达谱矩阵,然后取其平均值,得到整体细胞类型特异性基因表达谱矩阵 V V V。
网络结构
对于编码器网络 f Z ( ⋅ ) f_{Z}(\cdot) fZ(⋅),用于提取共享的隐空间中点的表示,STitch3D采用GAT将构建的3D空间graph融合。GAT包含一个图注意力层和一个dense层, f Z ( ⋅ ) f_{Z}(\cdot) fZ(⋅)图注意力层的输入是全局邻居graph A A A和归一化对数变换后的基因表达 X X X。 X X X的维数是被选中的高变基因的数量 G G G,GAT输出的维数为512。然后 f Z ( ⋅ ) f_Z(⋅) fZ(⋅)中的dense层以512维隐向量为输入,输出128维隐表示 Z Z Z。对于产生细胞类型比例的网络 f β ( ⋅ ) f_β(⋅) fβ(⋅),STitch3D采用单层dense网络将128维潜在表征 Z Z Z映射到 C C C维输出,其中 C C C为细胞类型的数量,并使用softmax激活。STitch3D采用 f α ( ⋅ , s ) f_α(⋅,s) fα(⋅,s)的单层dense网络,该网络以潜在表征 Z Z Z和切片标签为输入,生成slice-spot特异性效应 α n s α^s_n αns。对于正则化器 R A E R_{AE} RAE中的切片特异性解码器网络 f X ( ⋅ , s ) f_{X}(⋅,s) fX(⋅,s),使用维数为512的GAT网络。具体来说,图注意层的输入是全局邻接矩阵 A A A、潜在表示 Z Z Z和切片标签。输出到图注意层的维数为512。然后使用dense层,以512维隐藏向量为输入,重构基因表达矩阵。
为了展示如何利用STitch3D中的空间位置信息,作者描述了编码器网络 f Z ( ⋅ ) f_Z(⋅) fZ(⋅)和解码器网络 f X ( ⋅ , s ) f_X(⋅,s) fX(⋅,s)中使用的图注意层的细节。设spot i i i的嵌入 u i i n u_{i}^{in} uiin作为图注意层的输入。图注意层通过 A = [ A i , j ] ∈ R N × N A=[A_{i,j}]\in\mathbb{R}^{N\times N} A=[Ai,j]∈RN×N更新spot的embedding: u i o u t = σ ( ∑ j = 1 N a i , j W u j i n ) a i , j = A i , j e x p ( e i , j ) ∑ h = 1 N A i , h e x p ( e i , h ) e i , j = σ ( v 1 ⊤ W u i i n + v 2 ⊤ W u j i n ) u_{i}^{out}=\sigma(\sum_{j=1}^{N}a_{i,j}Wu_{j}^{in})\\a_{i,j}=\frac{A_{i,j}exp(e_{i,j})}{\sum_{h=1}^{N}A_{i,h}exp(e_{i,h})}\\e_{i,j}=\sigma(v_{1}^{\top}Wu_{i}^{in}+v_{2}^{\top}Wu_{j}^{in}) uiout=σ(j=1∑Nai,jWujin)ai,j=∑h=1NAi,hexp(ei,h)Ai,jexp(ei,j)ei,j=σ(v1⊤Wuiin+v2⊤Wujin)其中, W W W, v 1 v_{1} v1和 v 2 v_{2} v2为图注意层的网络参数。对于编码器网络 f Z ( ⋅ ) f_{Z}(\cdot) fZ(⋅), u i i n u_{i}^{in} uiin是基因表达级别的数据 X i X_{i} Xi, u i o u t u_{i}^{out} uiout为512维的hidden向量,对于解码网络 f X ( ⋅ , s ) f_{X}(\cdot,s) fX(⋅,s), u i i n u_{i}^{in} uiin为 Z i Z_{i} Zi。
模型训练
STitch3D在模型训练中使用Adam算法的一种变体Adamax进行随机优化。默认情况下,STitch3D的优化步数设置为20,000步,学习率为0.002,计算运行平均值的系数 β 1 = 0.9 β1 = 0.9 β1=0.9, β 2 = 0.999 β2 = 0.999 β2=0.999,权重衰减参数 λ = 0 λ = 0 λ=0。STitch3D的训练过程保证了大型数据集目标函数的收敛性,例如具有17086个spot的小鼠大脑数据集(Molecular atlas of the adult mouse brain)。
基因表达的去噪或插补
STitch3D的细胞型反卷积结果还可以对低质量基因的表达水平进行去噪,并在ST数据集中对未测量的基因表达进行补偿。STitch3D首先根据scRNA-seq reference计算出细胞类型特定的基因表达向量 V ′ = [ V c ′ ] ∈ R C × 1 V ' = [V'_{c}]\in\mathbb{R}^{C×1} V′=[Vc′]∈RC×1,其与获得高变基因的细胞类型特定基因表达谱矩阵 V V V相同。然后,使用STitch3D为每个spot估计的细胞类型比例 β i \beta_{i} βi,去噪或插补后的表达为 ∑ c = 1 C β i , c V c ′ \sum_{c=1}^{C}\beta_{i,c}V'_{c} ∑c=1Cβi,cVc′。
评估
在真实数据上评估
作者采用AUC来评估细胞型反卷积结果的可靠性。这里以人类DLPFC数据分析为例来说明如何计算AUC分数。DLPFC的空间研究为其生成的DLPFC Visium切片提供了层注释。单细胞数据集包含皮质层中的十种兴奋性神经元亚型。在单细胞研究中发现神经元亚型Ext_8_L5_6表达皮层5-6层标记基因,表明其定位于5-6层。基于这些信息,可以检查细胞型反卷积是否成功地通过估计的细胞类型分布识别出5-6层中Ext_8_L5_6亚型的定位。为此,根据对DLPFC Visium切片的手动标注,作者将5-6层中的spot标记为真,其他spot标记为假,作为ground truth。然后采用ROC分析,利用估计的Ext_8_L5_6亚型的细胞类型比例进行二元分类任务来预测spot标记,得到ROC曲线。AUC评分计算的是ROC曲线下的面积,可以用来指示估计的细胞类型比例是否分布在正确的区域。作者也对其他神经亚型进行了相同的ROC分析。为了进行公平的比较,在用于评估AUC分数之前,将每种方法的输出归一化,使每个spot的所有细胞类型比例之和等于1。
在模拟数据上评估
在模拟的ST数据集中,spot的细胞类型组成是已知的。为了进行公平的比较,将每种细胞型反卷积方法的输出归一化,使每个spot的所有细胞型比例之和等于1。作者使用了性能指标,包括Pearson相关系数(PCC)、结构相似性指数(SSIM)、均方根误差(RMSE)和Jensen-Shannon距离(JS),来定量比较每个spot估计的细胞类型比例与实际情况。
基准研究
作者首先使用人类背外侧前额叶皮层(DLPFC)数据集评估STitch3D的空间domain检测性能。原始研究中注释了6层(L1-L6)和白质(WM)。当应用于single slice时,STitch3D稳定地恢复了层结构。当应用于多切片分析时,STitch3D在各切片之间产生一致的结果,有利于三维重建(补充图1)。与层注释相比,STitch3D的结果具有相似的模式,表明其可靠性(图2a,b)。为了进行定量评估,将手工注释作为GT,用ARI评估。与single slice结果相比,STitch3D的多片结果获得了更高的ARI评分,这表明STitch3D能够跨切片借用信息(图2e)。将STitch3D与BayesSpace、SpaGCN和STAGATE等空间域识别方法进行了比较。与STitch3D不同的是,这些方法在皮质层的检测方面表现出较差的能力(图2e)。STitch3D的3D空间区域识别的优势在于它在共享潜在空间中有效地整合了切片(图2d)。潜在表征在所有切片中显示一致的模式,具有从WM和L6到L1的轨迹。
接下来,评估了STitch3D的细胞型反卷积性能。作者使用了相应的细胞类型reference,其中注释了10种兴奋性神经元亚型。根据ROC评估了细胞型反卷积的可靠性。ROC曲线下面积(AUC)值越高,表示性能越可靠。与单层结果相比,多层分析中STitch3D获得了更高的AUC评分(图2f)。在STitch3D的多层切片结果中,观察到1-4层的Ex_8_L5_6比例降低且噪声较小(补充图4)。STitch3D的多片结果可以清晰地表征神经元亚型的3D空间分布(图2c)。作者还将STitch3D与RCTD、Cell2location、Tangram、DestVI和CARD等细胞型反卷积方法进行了比较。与STitch3D相比,这些方法的AUC评分较低(图2f)。其中,Cell2location还允许进行多层分析。然而,从单片和多片实验中的AUC分数可以看出,它显示出有限的跨片借用信息的能力。
在基准研究之后,作者还对细胞型反卷积方法进行了全面比较。作者在三种不同的场景中模拟ST切片。由于已知模拟数据集中的真实细胞类型比例,可以使用总体准确性评分定量测量结果(整合PCC,SSIM,RMSE和JS)。
- 补充图1:基于人背外侧前额叶皮层(DLPFC)数据集的空间域检测方法比较。作者将所有比较方法的空域检测结果可视化,并对GT层进行标注。在单片和多片设置中均使用了STitch3D。
- 图2a:STitch3D对DLPFC数据的三维空间域检测结果
- 图2b:手工标注的三维可视化
- 图2c:STitch3D重建的Ex_8_L5_6神经元亚型的三维空间分布。比例值大于20%的点显示为紫色
- 图2d:用UMAP可视化了STitch3D在共享潜在空间中学习到的spot表示,通过切片索引、STitch3D分配的域标签和手动标注的标签着色
- 图2e:空间域检测方法的ARI比较
- 图2f:十种区域特异性神经元亚型的细胞型反卷积方法的AUC分数的小提琴图
- 图2g:在一项基于seqFISH+数据的模拟研究中,比较方法的总体精度得分(整合PCC,SSIM,RMSE和JS)
- 图2h:在基于STARmap数据的模拟研究中,比较方法的总体精度分数
- 图2i:MERFISH数据模拟研究中比较方法的总体精度分数
- 补充图4:使用DLPFC的细胞型反卷积方法的比较。从切片151676上的细胞型反卷积方法估计出的ex8_L5_6神经元亚型比例。
空间域检测方法本质还是聚类,但是这些聚类方法总是呈现出一个问题:在部分切片上聚类表现很好,但是换个切片,即使切片距离很近,聚类效果就呈现出很大差异,STitch3D不仅可以在不同切片的聚类中表现稳定,还能判断切片中每个spot的细胞类型占比。
实验
重建成年小鼠大脑
在本节中,作者展示了STitch3D可以准确地重建一个复杂的3D成年小鼠大脑。作者使用了横跨前后轴的35个冠状面切片和包含59个细胞类型的细胞类型reference。这里的3D重建任务是具有挑战性的,因为它需要在几十个切片中考虑批次效应,并区分细微的细胞亚型之间的差异。
- 补充图15:STitch3D对成年小鼠大脑的空间域检测结果。a为原始数据的UMAP图,用切片indices着色。b为由聚类标签着色的STitch3D学习表征的UMAP图和PAGA算法推断的空间轨迹。
- 图3a:在二维ST切片上可视化的STitch3D空间域检测结果中的簇1、2和5
- 图b和c:使用ICP和Allen Mouse Brain Common Coordinate Framework (CCFv3) 在STitch3D的对齐3D坐标中显示簇1,2和5
- 图d:在二维ST切片上显示的STitch3D细胞型反卷积结果中,估计了四种海马神经元类型的比例
- 图e和f:使用ICP和CCFv3的四种海马神经元类型在STitch3D对齐的3D坐标中所占比例的3D可视化,显示了比例值大于20%的点
- 图g:比较方法的AUC分数用六种区域细胞类型测量
- 图h和i:z-score中基因Tle4 (h)和Rell1 (i)的去噪和原始表达水平
- 图j:在6个二维切片上显示了皮层(簇1、2、5)间的伪时间空间轨迹
- 图k:包含虚拟矢状切片的平面可视化
- 图l:虚拟矢状面的空间域,以及四种海马神经元类型的比例
在对齐切片后,应用了STitch3D。STitch3D基于整合的spot表示将大脑划分为组织良好的3D域(补充图15b)。例如,标记为簇1、簇2和簇5的三个层状结构域形成了等皮质区(图3a-c)。这些簇在轨迹推断中表现出很强的连通性(补充图15b)。这些簇内的伪时间分析显示,所有切片的皮质生成模式一致(图3j)。簇3和9对应海马和丘脑区域,它们沿前后轴变化,表明STitch3D能够保留切片之间的生物异质性。通过这个例子,我们也验证了STitch3D的批次效应建模的有效性。当作者将其从STitch3D中移除时,切片的整合程度较差,导致不同切片之间的空间域检测不一致。
利用reference中的细粒度细胞类型特征,STitch3D揭示了3D细胞类型分布。例如,它准确地重建了四种海马神经元类型的分布(图3e和f)。这些分布与Allen参考图谱-小鼠大脑中标注的海马区CA1、CA2、CA3和DG正确匹配。为了进行定量比较,作者使用估计的四种海马神经元类型的比例来恢复CA1、CA2、CA3和DG区域,并使用ROC分析将结果与区域注释进行比较。对2-3层和5-6层的两种皮层神经元类型进行了类似的分析(图3g)。与其他方法相比,STitch3D的单片分析已经显示出总体上更好的AUC评分。通过借用切片间的信息,STitch3D在多切片分析中获得了更好的准确性。
ST研究中的一个问题是基因表达测量中的大量噪声。STitch3D通过将估计的细胞类型比例与reference的细胞类型特征相结合,实现基因去噪。作者利用基因Tle4和rel1证明了这一功能,它们在原始数据中具有噪声模式。应用STitch3D后,基因表达模式变得清晰(图3h和i)。
STitch3D的3D重建结果进一步提供了原始切片未捕获的组织的辅助视图。在使用冠状面切片重建的3D模型中,通过引入平行于前后轴的平面,创建了一个虚拟矢状面切片(图3k)。通过确定空间域和估计细胞类型比例,虚拟切片正确显示了矢状脑结构(图3l)。