引言

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

动态建模

我们采用广义动态模型来解析完整的转录动态过程。

通过这种方式，我们能够获得一些额外的发现，比如潜在时间以及潜在驱动基因的鉴定。

就像在之前的文章中一样，我们将这种方法应用于胰腺内分泌发育的研究。

# 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.pancreas()

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

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

动态模型

我们运行动态模型，以学习剪接动力学的完整转录动态。

该模型在基于似然的期望最大化框架中求解，通过迭代估计反应速率参数和细胞特异性潜在变量（即转录状态和细胞内潜在时间），从而学习每个基因的未剪接/剪接相轨迹。

scv.tl.recover_dynamics(adata)

scv.tl.velocity(adata, mode='dynamical')

scv.tl.velocity_graph(adata)

运行动态模型可能耗时较长。因此，您可以考虑将结果存储起来以便后续重复使用，使用 adata.write('data/pancreas.h5ad') 进行存储，稍后可以用 adata = scv.read('data/pancreas.h5ad') 来读取。

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

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

动力学速率参数

无需实验数据即可估算 RNA 转录、剪接和降解的速率。

这些参数有助于更好地理解细胞身份和表型异质性。

df = adata.var
df = df[(df['fit_likelihood'] > .1) & df['velocity_genes'] == True]

kwargs = dict(xscale='log', fontsize=16)
with scv.GridSpec(ncols=3) as pl:
    pl.hist(df['fit_alpha'], xlabel='transcription rate', **kwargs)
    pl.hist(df['fit_beta'] * df['fit_scaling'], xlabel='splicing rate', xticks=[.1, .4, 1], **kwargs)
    pl.hist(df['fit_gamma'], xlabel='degradation rate', xticks=[.1, .4, 1], **kwargs)

scv.get_df(adata, 'fit*', dropna=True).head()

基因特异性参数包括转录速率（fit_alpha）、剪接速率（fit_beta）、降解速率（fit_gamma）、切换时间点（fit_t_）、用于调整未剪接读数不足的缩放参数（fit_scaling）、未剪接和剪接读数的标准差（fit_std_u，fit_std_s）、基因似然度（fit_likelihood）、推断的稳态水平（fit_steady_u，fit_steady_s）及其相应的p值（fit_pval_steady_u，fit_pval_steady_s）、整体模型方差（fit_variance），以及用于将基因特异性的潜在时间对齐到一个通用的、基因共享的潜在时间的缩放因子（fit_alignment_scaling）。

潜在时间

动态模型可恢复潜在的细胞过程的潜在时间。这种潜在时间代表细胞的内部时钟，仅根据转录动态，近似细胞分化过程中经历的实际时间。

scv.tl.latent_time(adata)

scv.pl.scatter(adata, color='latent_time', color_map='gnuplot', size=80)

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

scv.pl.heatmap(adata, var_names=top_genes, sortby='latent_time', col_color='clusters', n_convolve=100)

Top-likelihood 基因

驱动基因表现出显著的动态行为，可通过在动态模型中以高似然度进行表征而被系统性地检测到。

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

scv.pl.scatter(adata, basis=top_genes[:15], ncols=5, frameon=False)

var_names = ['Actn4', 'Ppp3ca', 'Cpe', 'Nnat']

scv.pl.scatter(adata, var_names, frameon=False)

scv.pl.scatter(adata, x='latent_time', y=var_names, frameon=False)

聚类特异性Top-likelihood基因

此外，可为每个细胞聚类计算部分基因似然度，以实现聚类特异性的潜在驱动因子的鉴定。

scv.tl.rank_dynamical_genes(adata, groupby='clusters')

df = scv.get_df(adata, 'rank_dynamical_genes/names')

df.head(5)

for cluster in ['Ductal', 'Ngn3 high EP', 'Pre-endocrine', 'Beta']:
    scv.pl.scatter(adata, df[cluster][:5], ylabel=cluster, frameon=False)