4.基于python的scRNA-seq轨迹推断-RNA速率

参考：
[1] https://github.com/Starlitnightly/single_cell_tutorial
[2] https://github.com/theislab/single-cell-best-practices

背景

计算伪时间唯一可以利用的信息就是细胞表示之间的相似性，进一步，我们想挖掘更多信息来弥补现有缺陷。RNA速度估计利用额外的剪切和未剪切测量，目的是获得更准确的拟时序估计。

细胞转录组谱的变化是由一系列事件触发的：广义上说，DNA被转录以产生所谓的未剪切前体信使RNA（pre-mRNA）。未剪切的pre-mRNA包含与翻译相关的区域（外显子）以及非编码区域（内含子）。这些非编码区域被剪切掉，即被移除，形成剪切、成熟的mRNA。尽管scRNA-seq技术不能捕捉多个时间点的转录组，但它们确实包含了区分未剪切和剪切mRNA所需的信息。

我们用RNA速率来描述剪切RNA的变化：$$\frac{du}{dt}=\alpha -\beta u,\frac{ds}{dt}=\beta u-\gamma s$$其中，$\alpha$是转录速率，$\beta$是剪切速率，$\gamma$是剪切RNA的降解速率。$u,s$分别代表未剪切和剪切的RNA。

单细胞测量属于快照数据，因此不能准确根据时间点绘制图谱。相反，经典的 RNA 速度方法依赖于研究每个基因的未剪接和剪接 RNA 的细胞特异性元组 $(u,s)$。这些元组的集合形成了所谓的相图。假设转录、剪切和降解速率恒定，相图呈现杏仁形状。上方的弧对应于诱导阶段，下方的弧对应于抑制阶段。然而，由于现实世界的数据充满噪音，绘制未剪切与剪切counts的图并不能恢复预期的杏仁形状。 相反，需要首先对数据进行平滑处理。 传统上，该预处理步骤包括在细胞间相似性图中对每个细胞相对于其邻居的基因表达进行平均。

相图可以验证RNA速率的估计准确性。估计RNA速率的模型一般是稳态模型，EM模型，深度学习模型。

首先加载数据，我们使用胚胎期(E)15.5小鼠胰腺发育的scRNA-seq数据集，我们可以观察剪切和未剪切的比例（这部分需要测量得到），这是估计RNA速率必须的：

from tqdm import tqdm
from multiprocessing import Lock
tqdm.set_lock(Lock())  # manually set internal lock
tqdm.write("test")

import scanpy as sc
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import omicverse as ov
ov.plot_set()

import cellrank as cr
import scvelo as scv

adata = sc.read("./data/endocrinogenesis_day15.5.h5ad")
scv.pl.proportions(adata)
print(adata)
"""
AnnData object with n_obs × n_vars = 2531 × 27998
    obs: 'day', 'proliferation', 'G2M_score', 'S_score', 'phase', 'clusters_coarse', 'clusters', 'clusters_fine', 'louvain_Alpha', 'louvain_Beta', 'palantir_pseudotime'
    var: 'highly_variable_genes'
    uns: 'clusters_colors', 'clusters_fine_colors', 'day_colors', 'louvain_Alpha_colors', 'louvain_Beta_colors', 'neighbors', 'pca'
    obsm: 'X_pca', 'X_umap'
    layers: 'spliced', 'unspliced'
    obsp: 'connectivities', 'distances'
"""

一般adata.X = adata.layers['spliced']或者adata.X = adata.layers['spliced']+adata.layers['unspliced']。

下一步，我们进行预处理，并筛选出没有足够剪切/未剪切计数的基因，标准化并记录转换数据，限制在高度可变的顶部基因：

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

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

稳态模型

估计RNA速度时，假设基因是独立的，此外，还假设速率是恒定的，以及所有基因都有一个单一的共同的剪切速率，因此模型称为稳态模型，稳态本身位于相位图的右上角（诱导阶段）和原点（抑制阶段），稳态模型使用线性回归拟合来估计稳态比率。然后，RNA速率被定义为与这个拟合的残差。

稳态模型违反的假设是：所有基因具有相同的剪切速率，这往往会导致错误的结果。

在稳态模型下计算RNA速度，使用scVelo的velocity函数，并设置mode="deterministic"，我们还可以将速度投影到UMAP上：

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

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

EM模型

EM模型不再假设已达到稳态，也不假设基因具有共同的剪切速率。此外，所有数据点都用于推断完整的参数，以及剪切模型的基因和细胞特定的潜在时间。该算法使用期望最大化（EM）框架来估计参数。

EM方法如下：

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

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")

深度学习模型

前面都是基于统计的模型，现在也对应出现了基于深度学习的模型veloVI，由于环境问题，这里就不演示了。

获得伪时间

以EM模型为例，我们从中推断伪时间：

scv.tl.latent_time(adata)
adata.obs['latent_time_EMmodel']=adata.obs['latent_time'].copy()

adata.obs['velocity_pseudotime_EMmodel']=adata.obs['velocity_pseudotime'].copy()

在UMAP中可视化：

ov.pl.embedding(adata,basis='X_umap',
                   color=['latent_time_EMmodel'],
                   frameon='small',cmap='Reds')