12.1 前言

单细胞数据集允许以高分辨率研究生物过程，例如早期发育。例如，虽然分析的是单个细胞而不是整个组织，但细胞表型的变化无法随着时间的推移进行跟踪。这一事实源于单细胞测序方案的局限性。对细胞进行测序后，它会被破坏，因此无法在以后的时间点再次测量其定义特征。值得注意的是，实验技术不仅无法测量不同时间的一般细胞概况，而且无法测量这些变化发生的速度。可以使用轨迹推理(trajectory inference，TI)领域的工具来及时恢复沿发展景观的位置。然而，经典的TI方法不提供任何定向的动态信息。此外，这些算法传统上不考虑转录组读数和相似性之外的信息。

12.2 RNA velocity模型

细胞转录组谱的变化是由一系列事件触发的：广义上讲，DNA 被转录产生所谓的未剪接前体信使 RNA (pre-mRNA)。未剪接的前mRNA包含与翻译相关的区域（外显子）以及非编码区域（内含子）。这些非编码区被剪接，即去除，以形成剪接的成熟mRNA。虽然单细胞RNA测序 (scRNA-seq) 方案无法在多个时间点捕获转录组，但它们确实包含了分离未剪接和剪接mRNA读数所需的信息。

12.3 胰腺内分泌发生pancreatic endocrinogenesis中的RNA velocity推断

作为如何推断RNA velocity的实际示例，我们分析了胰腺的内分泌发育[Bastidas-Ponce et al., 2019]。在该系统中，前内分泌细胞（Ductal、Ngn3 low EP、Ngn3 high EP、Pre-endocrine）发育成四种内分泌细胞类型（Alpha、Beta、Delta、Epsilon）。在这里，我们使用scVelo[Bergen et al., 2020]来推断RNA velocity。

12.3.1 配置环境

from pathlib import Path

import scanpy as sc
import scvelo as scv

12.3.2 常规设置

scv.settings.set_figure_params("scvelo")
DATA_DIR = Path("../../data/")
DATA_DIR.mkdir(parents=True, exist_ok=True)

FILE_PATH = DATA_DIR / "pancreas.h5ad"

12.3.3 加载数据

为了使用scVelo估计RNA velocity，需要将未剪接和剪接计数存储在AnnData的层槽中。我们建议将整个计数（即未处理的数据）传递到scVelo管道。

adata = scv.datasets.pancreas(file_path=FILE_PATH)
adata

AnnData object with n_obs × n_vars = 3696 × 27998
    obs: 'clusters_coarse', 'clusters', 'S_score', 'G2M_score'
    var: 'highly_variable_genes'
    uns: 'clusters_coarse_colors', 'clusters_colors', 'day_colors', 'neighbors', 'pca'
    obsm: 'X_pca', 'X_umap'
    layers: 'spliced', 'unspliced'
    obsp: 'distances', 'connectivities'

12.4 数据预处理

由于scRNA-seq数据充满噪声且稀疏，因此必须对数据进行预处理，以便使用稳态或EM模型推断RNA velocity。第一步，我们过滤掉未剪接和剪接RNA均未充分表达的基因（此处至少20个）。接下来，对未剪接和剪接RNA的细胞大小进行归一化，并在adata.X log1p中进行计数转换，以减少异常值的影响。接下来，我们还识别和过滤高度可变的基因（此处为2000）。

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

Filtered out 20801 genes that are detected 20 counts (shared).
Normalized count data: X, spliced, unspliced.
Extracted 2000 highly variable genes.
Logarithmized X.

到目前为止，数据预处理与经典的scRNA-seq工作流程类似。对于RNA velocity，我们还通过其附近的平均表达来观察结果。 这可以使用scVelo的moments函数来完成。

sc.tl.pca(adata)
sc.pp.neighbors(adata)
scv.pp.moments(adata, n_pcs=None, n_neighbors=None)

computing moments based on connectivities
    finished (0:00:00) --> added 
    'Ms' and 'Mu', moments of un/spliced abundances (adata.layers)

在经典的工作流程中，我们会对数据进行聚类，推断细胞类型，并以二维嵌入方式可视化数据。幸运的是，对于胰腺数据，这些信息已经被先验计算并可以直接使用。

scv.pl.scatter(adata, basis="umap", color="clusters")

12.4.1 RNA velocity推断-稳态模型

第一步，我们计算稳态模型下的RNA velocity。在这种情况下，我们用mode="deterministic"来调用scVelo的velocity函数。

scv.tl.velocity(adata, mode="deterministic")

computing velocities
    finished (0:00:00) --> added 
    'velocity', velocity vectors for each individual cell (adata.layers)

虽然我们不鼓励过度解释高维速度向量到数据低维表示的投影，但scVelo提供了一种简单的方法。

scv.tl.velocity_graph(adata, n_jobs=8)
scv.pl.velocity_embedding_stream(adata, basis="umap", color="clusters")

computing velocity graph (using 8/96 cores)

    finished (0:00:05) --> added 
    'velocity_graph', sparse matrix with cosine correlations (adata.uns)
computing velocity embedding
    finished (0:00:00) --> added
    'velocity_umap', embedded velocity vectors (adata.obsm)

12.4.2 RNA velocity推断-EM模型

为了用EM模型计算RNA velocity，需要首先推断剪接动力学参数。scVelo的recover_dynamics函数负责完成。

scv.tl.recover_dynamics(adata, n_jobs=8)

recovering dynamics (using 8/96 cores)

    finished (0:01:00) --> added 
    'fit_pars', fitted parameters for splicing dynamics (adata.var)

拼接模型的参数是通过最大化给定的可能性来推断的。为了研究哪些基因最适合scVelo，我们可以研究相应的相图以及推断的轨迹（以紫色绘制）和稳态比率（紫色虚线）。此处，显示的五个基因中的三个（Pcsk2、Top2a、Ppp1r1a）呈现出（部分）杏仁形状的相图。我们观察到单一细胞类型（Top2a、Ppp1r1a）内或跨多种细胞类型（Pcsk2，从前内分泌细胞到 Alpha 和 Beta）的明显转变。就Nfib而言，我们观察到两个处于稳态的细胞群。这很可能是对Ngn3低/高EP细胞周围表型流形进行欠采样的人为因素。类似地，Ghrl在Epsilon细胞中高度表达，但由于簇大小较小，表达量很少。虽然当前的最佳实践仅限于手动分析模型拟合度和置信度，但最近提出的方法可以帮助自动化该过程（新方向）。 在这里，Nfib和Ghrl 将被分配较低的置信度分数。

top_genes = adata.var["fit_likelihood"].sort_values(ascending=False).index
scv.pl.scatter(adata, basis=top_genes[:5], color="clusters", frameon=False)

scv.tl.velocity(adata, mode="dynamical")
scv.tl.velocity_graph(adata, n_jobs=8)
scv.pl.velocity_embedding_stream(adata, basis="umap")

computing velocities
    finished (0:00:02) --> added 
    'velocity', velocity vectors for each individual cell (adata.layers)
computing velocity graph (using 8/96 cores)
    finished (0:00:02) --> added 
    'velocity_graph', sparse matrix with cosine correlations (adata.uns)
computing velocity embedding
    finished (0:00:00) --> added
    'velocity_umap', embedded velocity vectors (adata.obsm)

基于2D投影，EM模型更稳健地捕获导管细胞中的细胞周期。此外，稳态模型的投影表现出从阿尔法细胞到前内分泌细胞的“回流”。然而，为了进行严格的定量分析，我们建议使用CellRank等下游工具来评估模型差异并得出结论。