引言
本系列讲解 单细胞（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）、切换时间点（fitt）、用于调整未剪接读数不足的缩放参数（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)
