话不多说，开始我们的RNA velocyto

Introduction

single cell 技术最重要的优势：单细胞水平的分辨率 + 高通量带来的大样本

首先单细胞水平的分辨率使得 single cell 技术可以直接得到单细胞的表达谱，不需要像 bulk 时代那样需要取一堆细胞然后才能测得平均表达谱，因此在研究细胞分化上有天然的优势。并且可以对发育早期的少量细胞进行测序，得到异质性的结果。

v2-3b1d9c6f4944e6c02eddefadc2a5f75c_720w.jpg

高通量技术带来的大样本，使得可以更完整的刻画 cell state，并且对于各种需要进行统计推断的方法提供巨大的优势。

但是之前大多数 single cell 技术都是静态的刻画 cell state（即都是snapshot），没有能够引入细胞随时间的变化，living cell imaging tracking 的技术虽然能提供时间信息，但是能够覆盖到gene 数量是非常少的（通常也就不到10）。

而发育分化过程中 cell state 是动态变化的，因此怎么能够理解 cell state transition，怎么将这些 snapshot 拼成 video 就是一个很有意思且很重要的问题。

RNA velocity

18年 “RNA velocity of single cells” 首先引入 RNA velocity 这个概念，充分利用现有的数据来得到全新的概念。可以说是非常精彩了。

图片.png

核心的思想就是，利用 splicing 这个生物过程来引入时间这个要素，从而推断出 cell state 的 velocity（速率）。

具体是什么意思呢，首先从一个例子说起，当我们看到下图这些跑步中的中间状态的时候，我们脑海中是可以串联起这四个动作的。

图片.png

我们可以从其中一个 snapshot 知道图中小人的运动状态，并且预测运动的方向。本质原因是因为我们对跑步这个过程有一个动力学的认识。

那么我们回到细胞内的转录，我们只需要通过中心法则中依赖于时间的生物过程，并对这个过程进行动力学刻画即可以得到这些 snapshot 的未来方向，那么应该使用什么生物过程呢？首先需要满足两个条件，1. 依赖于时间，2. 现有的数据中应该可以直接得到

答案就是 splicing！ 首先 splicing 肯定依赖于时间，其次现有的测序 reads 天然的可以 mapping 到 intro 和 exon, 从而得到 unspliced 和 spliced 的RNA。

图片.png

Gene 在转录成 mRNA 的过程中会经过 splicing，而 splicing 的动力学我们是通过 ODE 来 model 的，因此我们便可以通过 unspliced mRNA (u)和 spliced RNA(s) ，以及动力模型( ODE model) 来计算 cell state 的 velocity 。我们可以得到

图片.png

上述便是 RNA velocity 的 key idea，至于随机模型，那就是考虑到真实值是离散的，所以利用化学主方程（Master equation）来研究，并且分析。

Beyond RNA velocity

• La Manno, Gioele, et al. "RNA velocity of single cells."Nature, 560.7719 (2018): 494-498.开篇之作，引入了 RNA velocity 的概念，但是在参数估计等方面有值得改进的地方
• Bergen, Volker, et al. "Generalizing RNA velocity to transient cell states through dynamical modeling."Nature biotechnology, 38.12 (2020): 1408-1414.Fabian 组的文章，在参数估计上有了比较大的改进，核心的改进有丢掉了 general constant 的假设，利用 EM 来估计 gene specific 的parameter，从而使得估计更加准确
• Li, Tiejun, et al. "On the Mathematics of RNA Velocity I: Theoretical Analysis."bioRxiv(2020).李铁军老师的文章，系统的梳理 RNA velocity 背后的理论基础，包括确定性模型(ODE)，随机模型(化学主方程)，连续化等拓展，并且之后其实还有分析不同模型产生的结果的文章，但 in preparation，还未能读到。
• Qiu X, Zhang Y, Yang D, et al. Mapping vector field of single cells[J]. Biorxiv, 2021非常精彩的文章，结合 metabolic labeling 的数据（可以标记新生成的 RNA）与 RNA velocity，从而对于 RNA velocity 有更真实意义上的刻画，并且将离散的 velocity 推广在cell state space 的 Vector Field，并且有一系列的分析.

上述几篇文章中可以称的上 beyond 的也就是 "Mapping vector field of single cells" (dynamo)，结合 metabolic labeling (能够标记新生成的 mRNA)，引入真实的时间，从而从两个时间尺度来刻画 velocity

包括三种model ，1. model I 单纯的 RNA velocity； 2. model II 单纯的 labeling RNA velocity，3. model III RNA velocity + labeling velocity。并且在离散的 velocity 基础上得到 vector field。

图片.png

接下来就是RNA velocity and beyond Theory 部分 ，主要介绍 model detail，怎么 estimate parameter等，最后写 beyond 的部分，也就是dynamo 的内容。

首先描述这个过程的模型（也就是如果知道那些 参数的情况下，mRNA 应该是怎么样的），然后说明怎么从已有的表达谱数据中估计这些参数。（本文只有最简单的利用 steady state 来估计） 这篇文章只讨论了确定性模型，没有从化学主方程的那种角度讨论，大概就相当于没有考虑这个离散性，不需要使用随机微分方程来刻画。

v2-39fb8de74af4e388c3e4355a407d3489_720w.png

RNA velocity 的关键在于利用 splicing 这一蕴含时间的动力学过程来刻画，从而在原来的 snapshot 中引入 velocity。

Deterministic Model

首先根据 RNA velocity 的 dynamics 可以得到相对应的 ODE

图片.png

对于转录本的生成过程我们可以知道，对于 ，我们有

图片.png

这个式子代表当转录启动的之后，转录的速率是一个常数 ，反之当转录关闭的时候，生成速率为 0。

而我们认为的 RNA velocity 也就是 spliced mRNA 的变化速度为

图片.png

在这个式子中其实已经做了 都是常数的假设，在这些假设的前提下，我们可以对上述方程求得显式解。

Case 1 不等于

解上述方程可以得到

图片.png

但是值得注意的是，上式代表的是一个 stage，如果 switch time 刚好被包含在这段计算的时间内，那么就需要分开处理，分阶段得到。 具体见下图也就是右边部分。

图片.png

分阶段

对于第一阶段 <= ，如果假设(₀,₀) = (0,0) ，我们可以用 上述的显式解写得

图片.png

而对于 第二个阶段也就是 < = 时，我们不妨定义(_s,_s) = ((),()) ，那么在接下来的 off 阶段我们有

图片.png

换句话就是两次利用方程的解析解，因为这个过程是被分成了两部分。

Case 2: 等于

和之前的区别大概就只是数值求解上的区别。所以对于 on stage 的解就是

图片.png

其中 <= ，如果此时 (₀,₀) = (0,0),可以得到这段时间的解为

图片.png

而在 off stage 可以得到

图片.png

而且其实也可以从之前的 不等于 的式子取极限得到。

Inference of RNA velocity

有了上面的模型，那么下一步就是通过数据来拟合这些参数。、、、。

从总体来说，我们可以知道、（,） (表示 switch time，switch 之后的状态)是无法从观测数据得到的，是 hidden variable，所以需要 EM 算法来求解。

Steady State Model

先从最简单的开始，如果我们认为 on stage 足够久的话，我们可以认为最后的状态是 steady state ，因此可以得到 steady state 的等式为

图片.png

而通过这个式子可以得到可以通过最小二乘法线性拟合来估计参数：

图片.png

最后一个等式假设 = 1（不失一般性的假设，也可以看做无量纲化，得到的相对值，这个假设也是开始的 Velocity 原文使用的) 据此可以得到 的估计。

图片.png

值得说明的是，在做回归估计的时候需要利用那些处在 steady state 的细胞（在 velocity 原文中也就是 那个 extreme 假设，只取 、图中两端的 cells）。

image

在图中的表现即只使用哪些绿色的点来估计斜率，可以看到紫色的估计会更接近真实值相比使用所有的细胞来做最小二乘线性拟合。

基于、 - *其实可以发现 velocity 的值其实就是、的一个线性组合，所以在这个情况下可能不能提供太多新的信息。

EM algorithm for the transient models

前面一小节似乎已经可以估计出参数了，但是为什么还要用 EM 算法呢？ 这是因为之前的 Steady state 的假设是不一定满足的，或者换言之，上一小节的估计需要找到那些已经处于 steady state 的细胞。（并不一定能在数据中找到，而且换言之就只利用了一部分数据来估计参数，那么自然就可能比较 sensitive）。

而利用 EM 算法来估计的原因是我们可以观测到的数据 (、)中无法知道一些hidden variable 、（，），对于每个 细胞和 每个 gene。 如果直接使用 MLE 之类的优化参数不现实。

因此采用 EM 算法迭代的优化。

大致的核心思想在于，首先随机初始化参数，然后借助参数可以得到对隐变量的估计，虽然再根据估计出的隐变量来重新估计参数，循环往复，左脚踩右脚逐渐上天hhh，但是需要说明的就是，如何重新估计参数是有讲究的，这一步保证了更新的 log-likelihood 是不减的。但也要注意的就是 EM 算法不保证能找到全局最优，因为 log likelihood 不减这个性质只能保证找到一个局部最优解。

生活很好，有你更好