RNA Velocity and Beyond 系列2—Model

分享是一种态度

时隔许久再次更新，本来想一口气写完模型部分还有参数估计部分，但是发现参数估计部分有点长，所以拆开写了。

本文主要参考李铁军老师挂在 biorxiv 上的文章:

Li, Tiejun, et al. "On the Mathematics of RNA Velocity I: Theoretical Analysis."bioRxiv(2020)

本文首先描述这个过程的模型（也就是如果知道那些参数的情况下，mRNA 应该是怎么样的），然后说明怎么从已有的表达谱数据中估计这些参数。（本文只有最简单的利用 steady state 来估计） 这篇文章只讨论了确定性模型，没有从化学主方程的那种角度讨论，大概就相当于没有考虑这个离散性，不需要使用随机微分方程来刻画。建议对这部分感兴趣的读者可以自行看李铁军老师的文章。没看过前文的可以考虑先看前文 RNA Velocity and Beyond 系列1—Introduction

Model

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

Deterministic Model

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

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

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

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

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

Case 1 β≠γ

解上述方程可以得到

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

分阶段

对于第一阶段 t≤ts，如果假设(u0,s0)≡(0,0)，我们可以用 上述的显式解写得

而对于 第二个阶段也就是 ts≤t 时，我们不妨定义(u0,s0):≡(u(ts),s(ts))，那么在接下来的 off 阶段我们有

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

Case 2:β＝γ

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

其中t≤ts，如果此时(u0,s0)≡(0,0)可以得到这段时间的解为

而在 off stage 可以得到

而且其实也可以从之前的 β≠γ 的式子取极限得到。

Inference of RNA velocity

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

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

Steady State Model

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

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

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

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

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

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

EM algorithm for the transient models

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

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

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

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

未完待续，具体 EM 估计参数的细节见下文