4.2 RNA velocity¶
Motivation¶
单细胞测序技术使我们能够以高分辨率研究生物过程，例如早期发育。然而，与整体组织分析不同，单细胞测序并不能追踪单个细胞表型特征随时间的变化。这是因为单细胞测序的实验过程会破坏细胞，一旦细胞被测序，它就会被破坏，因此无法在以后的时间点再次测量其特征。值得注意的是，现有的实验技术不仅无法测量不同时间点的整体细胞谱，也无法测量这些变化发生的速度。
- 轨迹推断 (TI) 领域可以帮助我们利用现有的数据估计细胞在发育过程中的位置。
然而，传统的 TI 方法无法提供定向的动态信息。此外，这些算法通常仅考虑转录组读段和相似性等信息，不包含其他方面的的信息。
Modeling RNA velocity¶
细胞转录组谱的变化是由一系列事件触发的：简单来说，DNA 会被转录成未剪切的前体信使 RNA (pre-mRNA)。未剪切的 pre-mRNA 包含翻译所需的区域 (外显子) 以及非编码区域 (内含子)。这些非编码区域会被剪切出 (即去除) 以形成剪切成熟的 mRNA。
虽然单细胞 RNA 测序 (scRNA-seq) 协议无法捕获多个时间点的转录组，但它们确实包含区分未剪切和剪切的 mRNA 读段的信息 {cite}`velo:LaManno2018, velo:Srivastava2019, velo:He2022, velo:Melsted2021}。
这一信息揭示了转录过程的动态性。例如，检测到更多的未剪切读段表明细胞正在积极转录，而检测到更多的剪切读段则表明细胞正在减少转录并可能更专注于蛋白质合成。
识别未剪切和剪切的读段允许我们构建一个描述剪接动力学的动力学模型 {cite}velo:Zeisel2011，并基于单细胞数据推断相应的模型权重。模型描述的剪切 RNA 变化称为 RNA 速度 {cite}velo:LaManno2018。当前的 RNA 速度模型假设基因特异的模型为：
$$ \begin{aligned} \frac{du_g}{dt} &= \alpha_g - \beta_g u_g\\ \frac{ds_g}{dt} &= \beta_g u_g - \gamma_g s_g, \end{aligned} $$
其中 $\alpha_g$ 表示转录速率，$\beta_g$ 表示剪接速率，$\gamma_g$ 表示剪切 RNA 的降解速率。虽然每个基因的动力学都是独立建模的，但为了简化记号，我们省略下标 $g$。尽管动力学系统参数估计领域已经研究得比较成熟，但推断算法需要知道每个观测数据关联的时间点。因此，这些传统方法不能用于推断 scRNA-seq 数据中的 RNA 速度及其模型参数。
Parameter inference¶
单细胞测序数据是快照数据，因此无法直接绘制成时间序列图。经典的 RNA 速度方法则采用另一种策略，研究每个基因的未剪切 $(u)$ 和剪切 $(s)$ RNA 读段组成的细胞特异性数据对 $(u, s)$。所有这些数据对的集合构成了所谓的相图 (phase portrait)。
假设转录、剪接和降解速率保持恒定，那么相图将呈现杏仁形状。上弧线对应于启动期，下弧线对应于抑制期。然而，由于真实世界的數據存在噪音，直接绘制未剪切读段和剪切读段的散点图并不能得到预期的杏仁形状。因此，数据需要经过平滑处理。传统的方法是在细胞相似性图中，对每个细胞的基因表达进行邻居平均 (averaging the gene expression of each cell over its neighbors)。
The steady-state model¶
第一种估计 RNA 速度的方法假设基因独立，并且底层动力学服从上述模型。此外，它还假设：（1）动力学已经达到平衡；（2）速率保持恒定；（3）所有基因具有相同的剪接速率。由于第一个假设，以下我们将此模型称为“稳态模型”。稳态值本身位于相图的右上角（启动期）及其原点（抑制期）。基于这些极端分位数，“稳态模型”通过线性回归拟合来估计稳态比率。然后，RNA 速度被定义为拟合的残差。
尽管“稳态模型”在某些系统中可以成功恢复发育方向，但其本身的模型假设限制了它的适用性。该模型容易被违背的两个假设是：所有基因具有相同的剪接速率以及实验过程中观察到了平衡状态。因此，在这些情况下使用该模型会得出错误的推断结果。此外，“稳态模型”仅考虑了数据的一部分，并且只推断稳态比率，而没有推断每个模型参数。
The EM model¶
为了克服“稳态模型”的限制，研究人员提出了几种扩展方法。目前最流行的一种是 scVelo 中实现的“EM 模型” {cite}velo:Bergen2020。“EM 模型”不再假设系统达到了平衡状态或所有基因共享相同的剪接速率。此外，它利用所有数据点来推断整个参数集以及每个基因和细胞特异的剪接模型的潜在时间。该算法使用期望最大化 (EM) 框架来估计参数。E 步中找到的未观测变量包括每个细胞的时间和状态（启动期、抑制期或稳态期）。所有其他模型参数都在 M 步中进行推断。
尽管“EM 模型”不再依赖“稳态模型”的关键假设，因此适用范围更广，但推断出的 RNA 速度仍可能违背已有的生物学知识 {cite}velo:Bergen2021, {cite}velo:Barile2021。 造成这种推理失败的原因主要有两个方面：
一方面，“EM 模型”仍然假设速率保持恒定。因此，只要这个假设不成立，例如在红细胞成熟过程中 {cite}velo:Barile2021，推断就会出现错误。另一方面，和它的前辈一样，该模型也依赖于相图。因此，只要基因的相图不符合预期的形状，该算法就本质上无法应用并会得出失败的推断。
胰腺内分泌发生中的 RNA 速度推断¶
为了展示如何推断 RNA 速度，我们以胰腺的内分泌发育为例进行分析 {cite}velo:BastidasPonce2019。在这个系统中，前内分泌细胞（导管细胞、Ngn3 低表达 EP、Ngn3 高表达 EP、前内分泌细胞）会分化为四种内分泌细胞类型（Α 细胞、β 细胞、δ 细胞、ε 细胞）。这里，我们使用 scVelo {cite}velo:Bergen2020 软件包来推断 RNA 速度。
Environment setup¶
from pathlib import Path
import scanpy as sc
import scvelo as scv
General settings¶
scv.settings.set_figure_params("scvelo")
DATA_DIR = Path("../../data/")
DATA_DIR.mkdir(parents=True, exist_ok=True)
FILE_PATH = DATA_DIR / "pancreas.h5ad"
Data loading¶
为了使用 scVelo 估计 RNA 速度，需要将未剪切和剪切的计数存储在 AnnData 的 layers 槽中。我们建议将完整计数（即未经处理的数据）传递给 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'
Data preprocessing¶
由于单细胞 RNA 测序数据通常嘈杂稀疏，为了使用 稳态模型 或 EM 模型 推断 RNA 速度，数据需要进行预处理。第一步，我们会滤除未剪切和剪切 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.
迄今为止的数据预处理步骤与经典的单细胞 RNA 测序工作流程类似。但在 RNA 速度分析中，我们还需要通过计算相邻细胞的平均表达值来平滑观测数据。这可以通过使用 scVelo 的 moments 函数来完成。
- 这里需要注意的是，相邻细胞指的是根据细胞间的相似性定义的邻居关系。scVelo 会利用这些邻居关系来对每个细胞的基因表达进行局部平滑。
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)
在典型的单细胞 RNA 测序分析工作流程中，我们会对数据进行聚类分析以识别细胞类型，并通过降维手段将数据可视化在二维嵌入空间中。幸运的是，对于胰腺数据而言，这些信息已经事先被计算出来，可以直接使用。
scv.pl.scatter(adata, basis="umap", color="clusters")
RNA velocity inference - Steady-state model¶
作为第一步，我们使用稳态模型计算 RNA 速度。这里，我们使用 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)
尽管我们不建议过度解读将高维 RNA 速度向量投影到低维数据空间的结果，但 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)
0%| | 0/3696 [00:00<?, ?cells/s]
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)
RNA velocity inference - EM model¶
为了使用 EM 模型 计算 RNA 速度，首先需要推断剪接动力学参数。这个步骤可以利用 scVelo 的 recover_dynamics 函数完成。
scv.tl.recover_dynamics(adata, n_jobs=8)
recovering dynamics (using 8/96 cores)
0%| | 0/835 [00:00<?, ?gene/s]
finished (0:01:00) --> added
'fit_pars', fitted parameters for splicing dynamics (adata.var)
通过最大化给定的似然值来推断拼接模型的参数。为了研究哪些基因最适合由 scVelo 拟合，我们可以研究相应的相图以及推断的轨迹（紫色曲线）和稳态比率（紫色虚线）。这里，五个显示的基因中的三个 (Pcsk2, Top2a, Ppp1r1a) 的相图呈 (部分) 杏仁形。我们观察到单个细胞类型 (Top2a, Ppp1r1a) 内部或跨越多个细胞类型 (Pcsk2，从 Pre-endocrine 到 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)
估计了动力学速率 (存储为 adata.obs 的列 fit_alpha, fit_beta, fit_gamma) 后，我们可以计算速度和向我们二维 UMAP 嵌入的投影。
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)
0%| | 0/3696 [00:00<?, ?cells/s]
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)
基于二维投影，EM 模型 更真实地捕捉到了 Ductal 细胞中的细胞周期。此外，稳态模型 的投影显示了从 Alpha 到 Pre-endocrine 细胞的“回流”。然而，为了进行严谨的定量分析，我们建议使用 CellRank 等下游工具 cite 来评估模型差异并得出结论。
Key takeaways¶
为了理解 RNA 速率分析是否适用于给定的数据集，我们需要注意以下几点：
- 推断 RNA 速率需要研究的发育过程时间尺度与 RNA 分子的半衰期相当。例如，在胰腺内分泌发生过程中可以满足这一要求 cite，但对于阿尔茨海默病或帕金森病等长期疾病则不适用。类似地，RNA 速率分析不适用于稳态系统，例如缺乏 (成熟) 细胞类型之间任何转型的外周血单核细胞。
- 仅当底层模型假设 (近似) 成立时，才能稳健可靠地推断 RNA 速率。为了检查这些假设，可以研究相图以验证它们是否呈现预期的杏仁形。如果一个基因包含多个明显的动力学，则应谨慎应用 RNA 速率分析，并且可能需要将数据子集到单个谱系。
- 经典上，高维 RNA 速率向量通过将其投影到数据的低维表示上进行可视化。这种用于验证假设的方法可能会错误和误导，因为投影的速率流高度依赖于 (1) 包含的基因数量和 (2) 选择的绘图参数。此外，投影质量会在低维嵌入的边界处下降 cite。
New directions¶
虽然 RNA 速率分析已经成功应用于许多系统，但一些模型限制仍然存在。 违反模型假设可能会导致错误的结果 cite，并且将高维速度向量投影到数据的低维表示上可能会产生误导。 为了克服这些陷阱，已经开发了一些工具。 例如，CellRank cite 利用推断出的速度场来推断细胞可能的未来状态。 由于该算法在数据的高维表示上运行，因此可以避免嵌入上误导性的速度流。 相比之下，最近的一份出版物试图提高低维嵌入的质量 cite。
为了缓和当前 RNA 速率推断的假设，已经提出了几种新方法 [cite]('velo:Qiao2021', 'velo:MarotLassauzaie2022', 'velo:Chen2022', 'velo:Riba2022', 'velo:Gu2022', 'velo:Gu2022-PLMR', 'velo:Gayoso2022')。 例如，这些方法不再假设恒定速率 cite，可以使用原始计数 cite，或者在变分推理框架中重新制定推理方法，使估计值与不确定性相关联 cite。 此外，为了帮助理解 RNA 速率分析是否可以推断单个基因或整个数据集，已经提出了不同的程序 cite。
References¶
{bibliography}
:filter: docname in docnames
:labelprefix: velo