多尺度动力学模型表征空间细胞稳态|Nat. Methods

作者，Evil Genius

今天参考的文章在Spatial transition tensor of single cells（Nat. Methods）

空间转移张量(STT)，一种通过多尺度动力学模型使用信使RNA剪接和空间转录组来表征空间细胞稳态的方法。通过学习四维转移张量和空间约束随机游走，重建细胞状态特异性动力学和空间状态转变。

关于空间的velocyto，大家可以参考文章单细胞空间数据分析之velocity

• 单细胞角度：RNA velocyto方法模拟信使RNA (mRNA)的动力学，将细胞的未来剪接状态投射到scRNA-seq数据上，以揭示细胞命运决定的方向性，并改进轨迹推断、低维嵌入和基因调控网络推断。

• 空间角度：空间转录组学测量单个细胞或一小群细胞spot并有额外空间信息，允许在空间中分析异质细胞状态。其中可以参考文章10X空间转录组速率分析（Velocyto）之SIRV。

目前空转使用RNA velocyto存在的问题

1、当前的模型无法捕捉复杂空间组织中的多谱系或多亚稳定状态，因为剪接和非剪接转录物水平可能由于非线性基因调控或多细胞信号传导而发生差异。

2、mRNA剪接的时间尺度在几分钟或几小时内，在此期间，当前的RNA速率模型收敛到一个全局平衡，然而，细胞状态的转变可能跨越几天到几周。

3、目前主要的RNA velocyto方法只关注剪接速率，而忽略了与基因调控密切相关的非剪切速率，这可以进一步提供关于“‘attraction force”研究特定细胞状态的信息。

文章中分析到的数据包括单细胞、HybISS spatial dataset、10X visium、华大高精度数据。

相对于我之前分享的方法，这个方法唯一的缺点就是把HE底片给整没了。

多尺度细胞attractor理论提供了一个工具来模拟不同时间尺度和分辨率的动力学，以及解释多稳定状态。在该理论中，基因表达的时间变化及其相互调节被建模为由一组微分方程组成的动力系统。稳定的细胞类型对应于基因调控轻微扰动下动力系统的多个局部稳定固定点(即多稳定状态)，细胞的表达状态被“trapped”，高度可塑性的过渡细胞被建模为系统的“saddle point”，使细胞可以通过一定的方向进行状态转换。但是这种方法不能直接用于空间转录组。

文章提出了一种空间过渡张量(STT)方法，利用未剪接和剪接的mRNA计数来重建空间转录组数据中的细胞attractor，从而可以量化空间attractor之间的过渡路径以及分析单个过渡细胞。与具有一个全局平衡的线性RNA速率模型不同，STT假设在联合未剪接(U) -剪接(S)计数空间中存在多个attractor，细胞在attractor之间进行转换。构建了一个跨越细胞、基因、剪接状态和attractor的四维过渡张量，并将每个attractor basin与吸引子特异性数量相关联。STT通过迭代改进张量估计和分解张量诱导和空间约束的细胞随机游走，连接了局部基因表达和剪接动态之间的尺度以及attractor之间的全局状态转换。此外，STT对与多稳定表达模式相关的基因进行了排序，并对具有相似STT特性的途径进行了分类。通过研究非空间和空间数据集，证明了STT在揭示不同时空尺度上发生的细胞的多稳定attractor和过渡特性方面的独特能力。

Overview of STT

STT的输入是S和U计数的单细胞基因表达矩阵，以及细胞注释，作为对它们所属细胞状态的初步判断。此外，空间转录组数据也需要每个细胞(或spot)的空间坐标。STT通过参数估计和动力学分解之间的迭代，构造了形状为RNC×2×K×NG的attractor速率张量，其中NC表示细胞数量，NG表示基因数量，K表示吸引子数量。其他基于张量的动力学量，包括attractor中细胞的“身份”，转移概率和转移路径，随后在该构造中获得。

STT使用以下基因表达和剪接动态的随机模型

其中Ui和Si是基因i的未剪接计数和已剪接计数。非线性函数fi(t, S1，…，SNG)模拟了其他基因如何调节基因i的产生速率。系统可以拥有多个fixed points或attractors，代表不同的细胞状态。参数βi表示mRNA剪接速率，γi表示剪接后的mRNA降解速率。独立的维纳过程项Wi.t和Zi.t表示基因表达中的噪声。与剪接动力学相比，这种随机性可以在更长的时间尺度上诱导噪声诱导的多稳定attractors之间的细胞状态转换.

当大多数细胞位于与不同细胞状态相对应的多个attractor basins 中时，一小部分细胞在 saddle points上进行过渡(细胞分布的自然假设)，未剪接的mRNA产生项可以扩展并近似为其线性扩展，从而引入attractor依赖的mRNA转录率。这种扩展允许对参数进行鲁棒估计，并初始化每个细胞的attractor方向速率的分配，称之为过渡张量。

通过构造inner-product velocity kernel ，张量提供了与连续随机微分方程(SDE)渐近一致的细胞随机游走描述。结合基因表达相似度的高斯核和细胞空间坐标，构建的细胞随机walk使每个attractor中的细胞具有一致的速率、转移方向和相似的基因表达。此外，构建的随机walk 倾向于细胞更有可能在物理空间中向其他空间相邻的细胞过渡。通过对随机游动进行粗粒化和attractor层次分解，得到不同attractor的cells’ membership functions。在张量模型构建和随机游走分解之间的每次迭代中，更新的membership functions通过引入attractor不确定性来改进式函数中的参数估计。然后在迭代过程中确定其动态与US空间中的attractor性质最一致的基因。包括一个监视器模块，具有正则化和早期停止策略，可以控制提高迭代的鲁棒性。最后，将描述attractor细节的张量流线以及描述长时间跃迁的粗粒度跃迁路径投影到低维动态流形上，以显示细胞状态的跃迁。

Overview of STT

STT在恢复多稳态细胞状态中的基准测试

首先应用STT分析了两个基于模拟多稳定系统的合成数据集。在双稳态切换电路中，STT中吸引子上的平均速率流线比RNA速率流线和其他方法更清晰地显示了两个吸引子的结构。虽然scVelo和UniTVelo计算的RNA速率流线倾向于偏离吸引子位置，但STT流线向吸引子收敛，从而提供了更可解释的切换开lanscope表示。此外，STT计算熵值来区分fixed point附近的稳定细胞和跨越saddle points的过渡细胞。正如流线跃迁张量的两个分量所示，只有将未拼接量和拼接量放在一起考虑时，才能同时揭示两个attractor basins。尽管剪接的张量与描述吸引子之间转换的标准RNA速率一致，但未剪接的张量自然地引入了一种“attraction force”，将细胞“拉”向每个吸引子的中心，而与cellDancer的流线相比，细胞被吸引到吸引子内的“末端”。未拼接计数为细胞状态的吸引子提供了STT中“吸引”水平的测量。为了进一步对STT的准确性进行基准测试，比较了STT未拼接或拼接张量分量与模型中真实速率的余弦相似度，发现STT在估计拼接和未拼接速率方面都名列前茅。此外，当对数据集进行下采样时，STT的性能显示出良好的鲁棒性。

接下来，分析了上皮-间充质转化(EMT)过程中的模拟基因调控回路，其中三个吸引子，分别表示为上皮(E)，间充质(M)和中间细胞状态(ICS)，可能在某些参数范围内共存。与scVelo计算的RNA速率相比，STT平均速率明显恢复了这三个吸引子。总的来说，STT能够重建单细胞基因表达数据集中复杂的多稳态细节。

Benchmarking of STT in simulation datasets of toggle-switch and EMT circuits

STT突出了中间细胞状态（ICSs）在命运决定中的作用

接下来分析了人肺A549细胞系EMT诱导实验中的scRNA-seq数据，包括TGFB1治疗后前7天收集的时间序列snapshot。STT识别出三个吸引子，分别是E、ICS和M，它们与数据采集的时间点顺序一致。此外，主要在诱导后8 h或1 d收集的ICS吸引子附近的细胞具有更高的熵值，这表明该状态比上皮状态(0天和8 h)和间充质状态(3天后)更具可塑性。这与提出的癌症中中间上皮和/或间质状态的表型可塑性是一致的。

利用转移向量预测上皮-间质landscope中吸引子连接的转移路径，发现从E到M的转移概率通量总是经过ICS。换句话说，上皮细胞经历EMT从未直接转换到间充质状态，而是首先获得中间性状。未拼接和拼接计数通常表现出吸引子的多重稳定性。多稳定性得分高的基因在不同吸引子内的未剪接数和剪接数表达水平不同，并在E-ICS-M过渡过程中呈现渐变变化。虽然在差异基因表达分析中，排名较高的多稳定基因(如ITGA11)并没有被显著检测到，但它们被发现在促进EMT转变和肿瘤进展方面很重要。虽然剪接动力学的张量流线显示了通过ICS从E到M的总体方向，这也与UniTVelo的结果一致，但STT和cellDancer预测的未剪接计数以及联合U-S空间的基因表达动态都表明，细胞在EMT期间被“吸引”到ICS盆地。这也与基于张量诱导多稳定核的CellRank吸收概率分析相一致。总的来说，张量分量和全局转换路径分析突出了ICS作为一个独特的吸引子盆地，在EMT期间充当枢纽状态。

此外，将STT应用于血液和胰腺发育数据集，发现其具有解决复杂状态转换的能力，其多稳定性张量核与CellRank分析一致。

Multistability of EMT in A549 cell lines with TGFB1 induction

STT识别空间吸引子和路径相似性

接下来，将STT应用于小鼠大脑发育的空间数据集。为了丰富未拼接和拼接计数以获得更好的张量估计，我们使用SIRV算法在E10和E11的40 μm处估算一个原始空间数据切片。与仅基于细胞相似性的聚类相比，STT识别的吸引子与原始出版物中不同细胞状态的空间位置和大脑区域注释一致:同一吸引子内的细胞往往具有相似的空间坐标，属于同一区域。此外，该方法对空间扩散核的权值、吸引子初始化、多稳定性基因滤波和吸引子个数具有鲁棒性。前脑和后脑吸引子的局部跃迁张量与UniTVelo分析结果一致

为了评估张量流线的生物学意义，进行了通路特异性分析，以评估与细胞状态转换和通路调节相关的功能。利用 Kyoto Encylopedia of Genes and Genomes(KEGG)知识库，基于多稳定基因的张量相关性计算各通路之间的相似度。事实上，特定路径张量表现出明显的吸引子动力学。基于张量相关的路径隐嵌入聚类揭示了发育过程中路径间空间状态转换的功能相似性。已知tgf - β和WNT通路在胚胎发生过程中表现出 cross-talk and cooperate，它们在潜在嵌入中来自不同的cluster，它们的张量流线方向相反，特别是在中脑和前脑吸引子中。大脑发育的另外两个重要通路，河马和甲状腺激素信号通路，也来自不同的通路张量cluster，在中脑和前脑区域显示相反的流线。总的来说，STT为细胞状态的空间组织和发育过程中调节状态转变的途径之间的关系提供了动态信息。

Transition tensor analysis of HybISS mouse brain spatial transcriptomics dataset

STT揭示了鸡心的空间吸引子和谱系

将STT应用于10X Visium技术测量的空间分辨鸡心数据。分析集中在数据集的最后一个时间点，即第14天，此时心脏的发育已经完成，心脏发生的完整事件和明确的空间边界已经完成。

利用SIRV输入的未拼接和拼接计数,STT识别出5个空间吸引子。其中，吸引子2与原研究中的“瓣膜”区域重合，主要由成纤维细胞组成。吸引子0主要由右心室区域的细胞组成。吸引子1主要位于“心房”区域，由红细胞组成。而其余的吸引子(3,4)分布在几个连接的区域，它们都包括心肌细胞注释表型的细胞。动态流形揭示了与不同细胞系相关的离散吸引子。吸引子1和2包含成纤维细胞和红细胞的空间定位谱系，在张量流线中都表现出“吸引力”。相比之下，吸引子3和4(都含有心肌细胞)内张量的流线表明它们在空间中的短暂性，并显示出向心房区域转移的趋势，这也可以在未拼接的组件之间的“吸引力”中观察到。这可以部分解释为心房中存在另一组肌细胞。总的来说，观察到的与空间区域或细胞类型注释的一致性表明STT有能力解剖空间分辨吸引子。

STT阐明了特定区域的空间吸引子和稳定性

接下来，分析了用bin size 60处理的高分辨率Stereo-seq小鼠成年冠状半脑数据，揭示了具有多种生物学功能的神经元细胞的复杂结构域。直接应用STT可以发现几个区域特异性的空间吸引子，这些空间吸引子与大脑区域的功能注释非常一致。张量的收敛流线表明，在皮层亚区(吸引子4)和纹状体背区(吸引子10)等区域，基因表达动态的多重稳定性得到了很好的维持。流线在丘脑区域(吸引子8)所有张量分量中向外流动，表明其具有较高的可塑性。基于其张量动力学的路径嵌入表明，先前已知的相互作用通路如cGMP-PKG和钙信号通路48具有相似的张量动力学。这也表明cGMP-PKG与催产素途径不同，其中流线表明主要差异发生在杏仁核区域。总体而言，结果表明STT可以通过吸引子分析发现空间区域并量化其稳定性，即使在成熟组织中也是如此。

Transition tensor analysis of Stereo-seq mouse coronal hemibrain spatial data

结论

量化和建模未剪接和剪接计数的相对丰度已成为从scRNA-seq数据集中剖析细胞状态转变的有效机制方法。为了连接基因表达、mRNA剪接和细胞谱系动态之间的不同时间尺度，并研究这些状态的潜在吸引子，研究团队开发了STT，用于（1）构建吸引子相关的过渡张量，（2）分析概率过渡路径和过渡细胞，以及（3）推断导致细胞状态多稳态性的基因。通过在（1）过渡张量模型的参数推断和（2）张量诱导的随机动力系统的多尺度分析之间进行迭代计算过程来实现这一目标。
与RNA速率模型相比，STT在揭示基因表达和剪接动态的潜在吸引子以及量化它们之间的过渡方面是独特的。通过假设多稳态性，STT对初始状态规格或隐藏时间校正具有鲁棒性。细胞成员函数在估计过渡张量时量化过渡细胞，自然地连接了下游多尺度动力学分析。

为了识别过渡细胞和推断过渡路径，STT利用计算的过渡张量，而不是直接使用RNA速率，如CellRank或Dynamo。研究团队发现，多稳态过渡张量在下游分析中与吸引子假设更兼容。STT在张量构建和动态解剖之间的迭代方案更好地确保了这种自洽性。然而，由于吸引子假设未考虑振荡动态，STT需要改进以捕捉具有强细胞周期效应的数据集的非平衡特征。

基于速率核的细胞随机walk是连接STT中张量推断和动态分解模块的关键，允许不同尺度上更好地连接动力学。理论分析表明，不同的速率核选择会导致形式为常微分方程或随机微分方程的不同连续极限。在STT中，使用内积核构建的细胞随机游走被证明与基因表达的随机化学Langevin模型一致，而采用余弦核正确恢复速率场的方向性，用于可视化吸引子内的局部流线。除了微分方程模型，未来将STT表述在RNA速率的化学主方程框架中可能会有深远的研究价值。

最后看看示例代码，完整的在https://github.com/cliffzhou92/STT/tree/release/example_notebooks，这里我就看一下10X Visium的数据。

import stt as st
import scanpy as sc
import anndata
import scvelo as scv
import numpy as np
import pandas as pd

data_dir = '../data/'
adata = sc.read_h5ad(data_dir+'mouse_brain.h5ad')
sc.pp.neighbors(adata, use_rep = 'xy_loc',key_added ='spatial')

sc.tl.leiden(adata,resolution = 0.3)


sc.pp.neighbors(adata)
sc.tl.leiden(adata,resolution = 0.3,key_added = 'exp_leiden')
sc.pl.scatter(adata, basis='xy_loc', color='exp_leiden')

U = adata.layers['unspliced']
S = adata.layers['spliced']
if 'toarray' in dir(U):
    U = U.toarray()
    S = S.toarray()
X_all = np.concatenate((U,S),axis = 1)
adata_aggr = anndata.AnnData(X=X_all)
sc.tl.pca(adata_aggr, svd_solver='arpack')
sc.pp.neighbors(adata_aggr)

sc.tl.leiden(adata_aggr,resolution = 0.15)
adata.obs['joint_leiden'] = adata_aggr.obs['leiden'].values

adata_aggr.obsm['xy_loc'] = adata.obsm['xy_loc']
sc.pl.embedding(adata_aggr, basis='xy_loc', color='leiden')

STT dynamical analysis

# Initialize attractor assignment by region annotation
adata.obs['attractor'] = adata.obs['Region']
adata_aggr = st.tl.dynamical_iteration(adata,n_states = 8, n_iter = 15, weight_connectivities = 0.5, n_neighbors = 50,thresh_ms_gene = 0.2, spa_weight =0.3)


sc.pl.embedding(adata, basis="xy_loc", color=["attractor"])

sc.pl.embedding(adata, basis="xy_loc", color=["Region"])

Pathway Analysis

st.tl.compute_pathway(adata,adata_aggr,'KEGG_2019_Mouse')


fig = st.pl.plot_pathway(adata,figsize = (10,8),size = 100,fontsize = 12)
for ax in fig.axes:
    ax.set_xlabel('Embedding 1', fontsize=20)  # Adjust font size as needed
    ax.set_ylabel('Embedding 2', fontsize=20)  # Adjust font size as needed
fig.show()

adata.obsm['X_xy_loc'] = adata.obsm['xy_loc']
adata_aggr.obsm['X_xy_loc']=adata.obsm['xy_loc']
adata_aggr.obsm['X_xy_loc_aggr']=adata.obsm['xy_loc']
adata.obsm['X_xy_loc_aggr']=adata.obsm['xy_loc']
st.pl.plot_tensor_pathway(adata,adata_aggr, pathway_name = 'Wnt signaling pathway',basis = 'xy_loc')

st.pl.plot_tensor_pathway(adata,adata_aggr, 'TGF-beta signaling pathway',basis = 'xy_loc')

st.pl.plot_tensor_pathway(adata,adata_aggr, 'Hippo signaling pathway',basis = 'xy_loc')

Sensitivity Analysis

import matplotlib.pyplot as plt
adata.obs['attractor'] = adata.obs['Region']
sc.pp.neighbors(adata, use_rep = 'xy_loc',key_added ='spatial')
nrows =1
weights = [0.1,0.3,0.5,0.7,0.9]
ncols = len(weights)

fig, axes = plt.subplots(nrows, ncols, figsize=(14, 4))

for i,sw in enumerate(weights):
#adata.obsm['X_umap'] = adata_aggr.obsm['X_umap']
    adata_aggr = st.tl.dynamical_iteration(adata,n_states =8, n_iter = 15, weight_connectivities = 0.5,n_components = 21, n_neighbors = 100,thresh_ms_gene = 0.2, use_spatial = True,spa_weight = sw, thresh_entropy = 0.1)    
    ax = axes[i]
    sc.pl.embedding(adata, basis="xy_loc", color="attractor",show = False, ax = ax)
    ax.set_title('spatial_weight='+str(sw))
plt.tight_layout()
plt.show()

Sensitivity on Initial Condition

adata.obs['spa_leiden']=adata.obs['leiden']

init_key = ['Region','spa_leiden','exp_leiden','joint_leiden']
fig, axes = plt.subplots(nrows, ncols, figsize=(12, 4))
for i,key in enumerate(init_key):
    adata.obs['attractor'] = adata.obs[key]
    adata_aggr = st.tl.dynamical_iteration(adata,n_states = 8, n_iter = 15, weight_connectivities = 0.5,n_components = 21, n_neighbors = 50,thresh_ms_gene = 0.2, spa_weight =0.3, thresh_entropy = 0.1)
    ax = axes[i]
    sc.pl.embedding(adata, basis="xy_loc", color="attractor",show = False, ax = ax)
    ax.set_title(key)
plt.tight_layout()
plt.show()

Sensitivity on MS score filter

adata.obs['attractor'] = adata.obs['Region']
adata_aggr = st.tl.dynamical_iteration(adata,n_states = 8, n_iter = 15, weight_connectivities = 0.5,n_components = 21, n_neighbors = 50,thresh_ms_gene = 0.2, spa_weight =0.3, thresh_entropy = 0.1)
sc.pl.embedding(adata, basis="xy_loc", color="attractor")

sc.pl.embedding(adata, basis="xy_loc", color="attractor")

adata.obsm['X_xy_loc'] = adata.obsm['xy_loc']
adata_aggr.obsm['X_xy_loc']=adata.obsm['xy_loc']
adata_aggr.obsm['X_xy_loc_aggr']=adata.obsm['xy_loc']
adata.obsm['X_xy_loc_aggr']=adata.obsm['xy_loc']
st.pl.plot_tensor(adata, adata_aggr, list_attractor = [3,5,6],basis = 'xy_loc',filter_cells = True, member_thresh = 0.1, density = 1)

Number of attractors

adata.obs['attractor'] = adata.obs['Region']
adata_aggr = st.tl.dynamical_iteration(adata,n_states = 8, n_iter = 15, weight_connectivities = 0.5,n_components = 21, n_neighbors = 50,thresh_ms_gene = 0.2, spa_weight =0.3, thresh_entropy = 0.1)
sc.pl.embedding(adata, basis="xy_loc", color="attractor")

adata_aggr.obsm['X_xy_loc']=adata.obsm['xy_loc']
scv.tl.velocity_graph(adata_aggr,n_jobs=-1,vkey = 'vj')
scv.pl.velocity_embedding_stream(adata_aggr, basis='xy_loc', color='attractor',vkey = 'vj',density = 2)

生活很好，有你更好