引言

本系列讲解 单细胞（scRNA-seq）中RNA“速率”分析教程[1]，持续更新，欢迎关注，转发！

差异动力学

在研究过程中，一个关键问题是我们需要处理包含多个谱系和过程的系统。在这种情况下，基因在不同亚群中往往会呈现出不同的动力学特征。由于不同的细胞状态和谱系通常受到基因调控网络不同变化的影响，因此它们可能会表现出不同的剪接动力学。这就使得一些基因在相空间中呈现出多样的轨迹。

为了解决这个问题，我们可以利用动态模型来进行差异动力学的似然比检验。通过这种方法，我们能够识别出那些其动力学行为无法通过一个单一的整体动态模型来很好地解释的细胞簇。随后，将细胞类型按照它们各自不同的动力学模式进行聚类，就可以分别对每种模式进行拟合。

为了更好地说明这一点，我们将差异动力学分析应用于齿状回神经发生的研究。齿状回神经发生过程包含了多个具有异质性的亚群。

# update to the latest version, if not done yet.
!pip install scvelo --upgrade --quiet

import scvelo as scv
scv.logging.print_version()

scv.settings.verbosity = 3  # show errors(0), warnings(1), info(2), hints(3)

scv.settings.presenter_view = True  # set max width size for presenter view

scv.settings.set_figure_params('scvelo')  # for beautified visualization

准备数据

数据处理主要包括基因选择、对数标准化以及计算矩。

adata = scv.datasets.dentategyrus()

scv.pp.filter_and_normalize(adata, min_shared_counts=30, n_top_genes=2000)

scv.pp.moments(adata, n_pcs=30, n_neighbors=30)

基础速度估计

scv.tl.velocity(adata)

scv.tl.velocity_graph(adata)

scv.pl.velocity_embedding_stream(adata, basis='umap')

差异动力学检验

不同的细胞类型和谱系可能会有不同的动力学特征，因为它们可能受到不同网络结构的调控。即使细胞类型或谱系之间存在联系，由于选择性剪接、选择性聚腺苷酸化以及降解过程的调节，动力学也可能存在差异。

动态模型可以帮助我们通过差异动力学的似然比检验来解决这一问题，从而识别出那些其动力学行为无法通过一个单一的整体动态模型来充分解释的细胞簇或谱系。每种细胞类型都会被检验，看是否通过独立拟合能够获得显著改善的似然度。

似然比遵循渐近卡方分布，可以进行显著性检验。需要注意的是，为了提高效率，默认情况下使用正交回归而不是完整的相轨迹来检验一个细胞簇是否能够被整体动力学很好地解释，或者是否表现出不同的动力学。

var_names = ['Tmsb10', 'Fam155a', 'Hn1', 'Rpl6']

scv.tl.differential_kinetic_test(adata, var_names=var_names, groupby='clusters')

scv.get_df(adata[:, var_names], ['fit_diff_kinetics', 'fit_pval_kinetics'], precision=2)

kwargs = dict(linewidth=2, add_linfit=True, frameon=False)

scv.pl.scatter(adata, basis=var_names, add_outline='fit_diff_kinetics', **kwargs)

以 Tmsb10 为例，内皮细胞表现出一种动力学行为（用黑色线条拟合并标出），这种行为无法通过整体动态（紫色曲线）很好地解释。

diff_clusters=list(adata[:, var_names].var['fit_diff_kinetics'])

scv.pl.scatter(adata, legend_loc='right', size=60, title='diff kinetics', add_outline=diff_clusters, outline_width=(.8, .2))

测试高似然基因

在筛选高似然基因的过程中，我们发现一些基因呈现出多种动力学模式。

scv.tl.recover_dynamics(adata)

#adata.write('data/pancreas.h5ad', compression='gzip')
#adata = scv.read('data/pancreas.h5ad')

top_genes = adata.var['fit_likelihood'].sort_values(ascending=False).index[:100]

scv.tl.differential_kinetic_test(adata, var_names=top_genes, groupby='clusters')

尤其是那些与主要颗粒细胞不同的细胞类型，例如 Cck/Tox、GABA、内皮细胞和小胶质细胞，它们频繁出现。

scv.pl.scatter(adata, basis=top_genes[:15], ncols=5, add_outline='fit_diff_kinetics', **kwargs)

scv.pl.scatter(adata, basis=top_genes[15:30], ncols=5, add_outline='fit_diff_kinetics', **kwargs)

重新计算速度

最后，可以利用多种竞争性动力学模式的信息来重新计算速度。

scv.tl.velocity(adata, diff_kinetics=True)

scv.tl.velocity_graph(adata)

scv.pl.velocity_embedding(adata, dpi=120, arrow_size=2, arrow_length=2)

[1]Source: https://scvelo.readthedocs.io/en/stable/DifferentialKinetics.html