实用指南：【论文阅读笔记】VeloCycle

文章目录

• • 流形
  • 似然
  • 贝叶斯模型
  • 变分分布（SVI）
  • 变分分布（LRMN）
  • 模型实现
  • 估算恒定细胞周期速度的近似点
  • 数据集

流形

• 利用单细胞RNA测序(scRNA-seq)，得到一个细胞的基因数$n\sim 10^4$

• $u_c$，$s_c$表示未剪切mRNA和已剪切mRNA的表达量，记$y_c=(u_c,s_c)$，$u_c,s_c\in {\mathbb{N}}^n$

• 每个细胞的潜在坐标$x$映射到凭借确定性函数$s(x)$（其中 s 表示“剪接”）描述的剪接基因表达水平流形$M$

• 通过根据障碍的生物学结构选择流形拓扑。例如，给定一个周期性过程，如细胞周期，我们能够取$x\in S_1$来体现周期性。

• 每个细胞$c$的测量值将通过真实的噪声模型与$M$上的相应位置相关联。在scRNA-seq的情形下，相关的噪声模型由负二项分布（NB）组成：
  $$\begin{gathered} Y_{gc}\sim NB[y_g(x_c),{\alpha }_g] \\ {y}_g(x_c)=\mathbb{E}[Y_{gc}]=(s_g(x_c),u_g(x_c)) \\ {\alpha }_g=({\alpha }_g^s,{\alpha }_g^u) \end{gathered}$$
  $c$是细胞，$g$是基因。这里假设了${\alpha }_g$与$x$独立。

• 在高维基因表达空间中，我们希望有一个描述RNA速度$\frac{d\tilde{s}}{dt}$的速率方程，该方程取决于剪接和非剪接RNA计数的预测：
  $$\frac{d{\tilde{s}}_g}{dt}=F({\tilde{s}}_g,{\tilde{u}}_g)={\beta }_g{\tilde{u}}_g-{\gamma }_g{\tilde{s}}_g$$
  其中，${\tilde{s}}_g$ 和 ${\tilde{u}}_g$分别是时间依赖的估计的剪接和未剪接RNA水平，${\beta }_g$ 和 ${\gamma }_g$是基因依赖的RNA剪接和降解速率。
  方程中的 $F$不显式依赖时间，剪接和降解速率被视为常数。

• 假设：存在一个自洽且确定性的方程描述潜在空间 x(t) 的动力学：
  $$\frac{dx}{dt}=V(x)$$
  其中 $V(x)$是潜在空间中的向量场。于是$\tilde{s},\tilde{u}$通过$x$传递成为时间依赖的：
  $$\begin{gathered} \tilde{s}(t)=s(x(t)) \\ \tilde{u}(t)=u(x(t)) \end{gathered}$$

• 由上述假设可能得到流形限制条件下的RNA速度：
  $$\frac{d{s}_g(x(t))}{dt}=({\nabla }_x{s}_g)\cdot V(x(t))={\beta }_g{u}_g(x(t))-{\gamma }_g{s}_g(x(t)),\forall g$$
  这里使用了链式法则。这个式子将左侧的低维流形拓扑与右侧的生物学连接起来。
  ${\beta }_g$和${\gamma }_g$是基因特异性剪接和降解率。值得注意的是，控制基因动力学的参数（$\beta$和$\gamma$）原则上也可能取决于$x$。
  ${\nabla }_{x}s$形成了切空间的$m$ 维基，$V(x(t))$提供了速度向量在该基中的分量。

• 我们可以依据该式子估算生物过程的实际持续时间：
  $$\Delta t_{s_0,s_1}={\int }_{{\Gamma }_{s_0}^{s_1}}\frac{1}{s}ds={\int }_{{\Gamma }_{x_0}^{x_1}}\frac{1}{V(x)}dx=\Delta t_{x_0,x_1}$$
  其中 ${\Gamma }_{x_0}^{x_1}$是连接两个点$x_0$ 和 $x_1$ 的轨迹 $x(t)$，并使用了轨迹变量$s(x)$ 的变化。

• 假设$M$在拓扑上是一个圆，将坐标$x$写成$\phi \in S^1$，于是动力学方程变成：
  $$\begin{gathered} \frac{d}{dt}{s}_g(\phi (t))=\frac{d}{d\phi }{s}_g(\phi )\omega (\phi )={\beta }_g{u}_g-{\gamma }_g{s}_g \\ E[s_{gc}]=s_g({\phi }_c)=\exp \left(\sum _f{v}_{gf}{\tilde{\zeta }}_f({\phi }_c)\right) \end{gathered}$$
  其中我们假设${\beta }_g$ 和 ${\gamma }_g$在细胞周期内是常数。

• 典型的细胞周期基因表现出只能用少数谐波描述的特征，因此，在展开中，大家将考虑$k$个傅里叶分量（在实践中，默认使用一个谐波）。又由于$s(\phi )$是正的，记：
  $$\begin{gathered} \log (s_g({\phi }_c))=\sum _f{v}_{gf}{\tilde{\zeta }}_f({\phi }_c) \\ {v}_g={\left(\begin{array}{cccccc}{a}_g^0& a_g^1& b_g^1& \cdots & a_g^k& b_g^k\end{array}\right)}^{\top } \\ \tilde{\zeta }(\phi )={\left(\begin{array}{cccccc}1& \cos (\phi )& \sin (\phi )& \cdots & \cos (k\phi )& \sin (k\phi )\end{array}\right)}^{\top } \end{gathered}$$
  这里$v_g$是用实数表示的基因傅里叶参数的向量。使用链式法则后得到$u(\phi )$：
  $$\begin{gathered} \frac{d}{dt}{s}_g(\phi (t))=\omega (\phi )s_g(\phi )\sum _f{v}_{gf}\frac{d}{d\phi }{\tilde{\zeta }}_f(\phi ) \\ \log (u_g(\phi ))=-\log ({\beta }_g)+\log (\omega (\phi )\sum _f{v}_{gf}{\partial }_{\phi }{\tilde{\zeta }}_f(\phi )+{\gamma }_g)+\log (s_g(\phi ))\forall g \\ E[u_{gc}]=u_g(\phi )=\frac{s_g(\phi )}{{\beta }_g}(\omega (\phi )\sum _f{v}_{gf}{\partial }_{\phi }{\tilde{\zeta }}_f(\phi )+{\gamma }_g) \end{gathered}$$

似然

• 可以根据剪接 RNA ($S_c$) 和未剪接 RNA ($U_c$) 的计数数据，计算每个细胞的似然函数。
• 全联合似然函数$P(\left\{(S_c,U_c)\right\}|\theta )$由以下表达式组成：
  $$\begin{gathered} P(\left\{(S_c,U_c)\right\}|\theta )=\prod _{gc}P(S_{gc},U_{gc}|\omega (\phi ),{\phi }_c,v_g,{\beta }_g,{\gamma }_g,{\alpha }_g) \\ P(S_{gc},U_{gc}|\theta )=P_s(S_{gc}|v_g,{\alpha }_g^s,{\phi }_c)\times P_u(U_{gc}|\omega (\phi ),{\beta }_g,{\gamma }_g,v_g,{\phi }_c,{\alpha }_g^u) \\ {P}_s(S_{gc}|\dots )=\text{NB}(s_g({\phi }_c)=F[v_g,{\phi }_c],{\alpha }_g^s) \\ P(U_{gc}|\dots )=\text{NB}(u_g({\phi }_c)=G[\omega ({\phi }_c),{\beta }_g,{\gamma }_g,v_g,{\phi }_c],{\alpha }_g^u) \end{gathered}$$
  其中$\theta$表示参数，$F,G$表示$s_g,u_g$与其他量的依赖关系。

贝叶斯模型

• 通过结合生物学定义的先验（priors）和从数据中得出的经验贝叶斯先验，近似计算联合后验概率分布$P(\theta |S_c,U_c)$：
  $$P(\theta |S_c,U_c)=\frac{P(S_c,U_c|\theta )P(\theta )}{P(S_c,U_c)}=\frac{{\prod }_{gc}P(S_{gc}|\theta )P(U_{gc}|\theta )P(\theta )}{\int {\prod }_{gc}P(S_{gc}|\theta )P(U_{gc}|\theta )P(\theta )d\theta }$$
  其中先验$P(\theta )$为：
  $$\begin{gathered} v{\omega }_t\sim \mathcal{N}([0,0,0],[3^2,0.05^2,0.05^2]) \\ \log ({\gamma }_g)\sim \mathcal{N}(0,0.5^2) \\ \log ({\beta }_g)\sim \mathcal{N}(2,3^2) \\ {\alpha }_g\sim \text{Gamma}(1.0,2.0) \\ {v}_{gt}\sim \mathcal{N}({\mu }_{gt}^v,{{\sigma }_{gt}^v}^2) \\ \phi x{y}_c=\text{ProjNormal}(\phi x_c,\phi y_c) \end{gathered}$$
• 通过经验贝叶斯（Empirical Bayes）设置以下参数：
  $$\begin{gathered} {\mu }_{gt}^v=[\log (\text{mean}c(Sgc)),0,0] \\ {\sigma }_{gt}^v=\left[\begin{array}{c}\frac{1}{2}\cdot {\text{std}}_c(Sgc+1)\\ \frac{1}{4}\cdot {\text{std}}_c(Sgc+1)\\ \frac{1}{4}\cdot {\text{std}}_c(Sgc+1)\end{array}\right] \\ \phi x_c=ϵ\cos ({\Phi }_c) \\ \phi y_c=ϵ\sin ({\Phi }_c) \end{gathered}$$
  其中${\Phi }_c={\tan }^{-1}({\omega }_{2c},{\omega }_{1c})$

变分分布（SVI）

• 变分分布 $P(\{v{\omega }_t,\{{\phi }_c\},\{v_{gt}\},\{{\beta }_g\},\{{\gamma }_g\},\{{\alpha }_g\})$被分解为多个独立分量的乘积，其形式为：
  $$P(\{v{\omega }_t,\{{\phi }_c\},\{v_{gt}\},\{{\beta }_g\},\{{\gamma }_g\},\{{\alpha }_g\})=\prod _{c}P(v{\omega }_t)P({\phi }_c)P(v_{gt})P({\beta }_g)P({\gamma }_g)P({\alpha }_g)$$
• 变分分布参数化如下：
  $$\begin{gathered} P(v{\omega }_t)\sim \mathcal{N}(\widehat{\mu v{\omega }_t},{\widehat{\sigma v{\omega }_t}}^2) \\ P(v_{gt})\sim \mathcal{N}(\widehat{{\mu }_{v_g}^v},\widehat{{\sigma }_{v_g}^v{}^2}) \\ P({\alpha }_g)=\text{Delta}(\widehat{{\alpha }_g}) \\ P(\log ({\gamma }_g))\sim \mathcal{N}(\widehat{{\mu }_{\log {\gamma }_g}},\widehat{{\sigma }_{\log {\gamma }_g}^2}) \\ P(\log ({\beta }_g))\sim \mathcal{N}(\widehat{{\mu }_{\log {\beta }_g}},\widehat{{\sigma }_{\log {\beta }_g}^2}) \\ P(\phi x{y}_c)\sim \mathcal{N}([\widehat{\phi x_c},\widehat{\phi y_c}],[1,1]) \end{gathered}$$
  

变分分布（LRMN）

• 低秩多变量正态（LRMN）模型考虑了观测数据之间的相关结构，基于变分推断（VI）构造的变分分布，观察到的联合后验由 MCMC 采样估计。具体而言，大家允许协方差和建立速度场$v_{\omega t}$以及动力学参数${\beta }_g$ 和 ${\gamma }_g$之间的关系。
• 后验因子分解如下：
  $$\begin{array}{rl}P\left(\{\nu {\omega }_t\},\{{\phi }_c\},\{{\nu }_{gt}\},\{{\beta }_g\},\{{\gamma }_g\},\{{\alpha }_g\}\right)& =P\left(\{{\gamma }_g\},\{\nu {\omega }_t\}\right)\prod _{g}P({\beta }_g\mid {\gamma }_g)P({\alpha }_g)\prod _{t}P(\nu {\omega }_t)P({\nu }_{gt})\prod _{c}P({\phi }_c)\end{array}$$
• ：就是具体公式
  $$\begin{gathered} \mathbf{x}\equiv \left[\log ({\gamma }_1),\log ({\gamma }_2),\dots ,\log ({\gamma }_{n_g}),\nu {\omega }_0,\nu {\omega }_1,\dots ,\nu {\omega }_{n_t-1}\right] \\ \mathbf{\Sigma }=\hat{\mathbf{F}}{\hat{\mathbf{F}}}^{\top }+\mathrm{diag}(\hat{\mathbf{d}})\text{where}\hat{\mathbf{F}}\in {\mathbb{R}}^{(n_g+n_t)\times k},\text{with}k=5 \\ P(\{\log ({\gamma }_g)\},\{\nu {\omega }_t\})=P(\mathbf{x})=\text{MultivariateNormal}(\hat{\mathbf{m}},\mathbf{\Sigma }) \\ {\mu }_{\log {\beta }_g|\gamma }={\hat{\mu }}_{\log {\beta }_g}+{\hat{\rho }}_g\cdot {\hat{\mu }}_{\log {\beta }_g}\cdot \frac{\log ({\gamma }_g)-{\hat{\mu }}_{\log {\gamma }_g}}{{\sigma }_{\log {\gamma }_g}}\text{with}{\hat{\rho }}_g\in [0,1] \\ {\sigma }_{\log {\beta }_g|\gamma }={\hat{\mu }}_{\log {\beta }_g}\sqrt{1-{\hat{\rho }}_g^2} \\ P(\log ({\beta }_g)\mid \log ({\gamma }_g))=\mathcal{N}({\mu }_{\log {\beta }_g|\gamma },{\sigma }_{\log {\beta }_g|\gamma }^2) \\ P({\phi }_c)=\mathcal{N}([\hat{\phi }{\mathbf{x}}_c,\hat{\phi }{y}_c],[1,1]) \\ P({\nu }_{gt})=\mathcal{N}({\hat{\mu }}_{gt}^{\nu },{\hat{\sigma }}_{gt}^{\nu 2}) \\ P({\alpha }_g)=\text{Delta}(\widehat{{\alpha }_g}) \end{gathered}$$

模型实现

• 模型实现旨在估算联合后验概率分布的近似值，涉及角细胞周期速度 ($v{\omega }_t$)，和 $S^1$流形上的参数（）。该实现分两个步骤进行：流形学习和速度学习。
• 在流形学习中，我们估计每个细胞沿细胞周期流形 ($\varphi$) 的位置，以及每个基因的傅里叶级数（$v$）。
• 所有变量初始化为先验的均值。先验均值经过以下两种方式确定：
• • 应用数据的前两个主成分 ($\varphi$)，这是一种降维方法，提取数据的低维结构。
• • 使用每个基因剪接表达量 ($v$) 的均值和标准差 ($s.d.$)，以反映基因表达的统计特性。
    剪接计数 (ElogS) 的期望值从真实素材和负二项分布 (NB) 建模得出，允许捕捉表达数据的离散性和过分散性。
• 为适应不同数据集或批次间平均表达水平的差异，模型引入了第一个基因谐波系数的偏移项 ($\Delta v$)。
• 速度学习的目标是基于流形学习的结果，估算傅里叶系数、角速度 ($v_{\omega }$) 以及速度动力学参数 ($\gamma$ 和 $\beta$)。
• 所有变量初始化为先验的均值。特别地：
• • 角速度 ($v\omega$) 的先验均值假设为零，反映了对零细胞周期速度的假设。
• • 其他变量（如傅里叶系数和动力学参数）也初始化为先验均值，具体取决于流形学习阶段的估计结果。
• • 为了确保模型输出满足生物学意义上的正值约束，独特是在方程 (10) 中 ($\omega (\varphi ){\sum }_f{v}_{gf}{\partial }_{\varphi }{\zeta }_f(\varphi )+{\gamma }_g$)，学习过程中引入了 ReLU 函数。
• 我们使用SVI求解VeloCycle模型，并应用ClippedAdam优化器和ELBO损失函数，从第一次到最终一次训练迭代，学习率从0.03衰减到0.005。
• 提供了提前终止选项：如果前 100 次迭代的均值损失与前 10 次迭代的均值损失相差小于五个单位，则停止训练。
• 速度动力学参数$\gamma$ 和 $\beta$受到生物学约束的限制：
• • ${\gamma }_g$的范围为 [0.5, 1.5] h${}^{-1}$，表示基因特定的动力学速率。
• • 周期 $T=2\pi /{\omega }_0$的范围为 [6, 50] h，反映细胞周期的生物学合理时间范围。
• 速度谐波系数的先验均值设为0，标准差为 3.0，反映了对无初始速度的假设，同时允许较大的变异性以适应内容变化。所有先验可以通过 ‘velocycle.preprocessing’ 套件中的函数修改，并通过 Pyro 模型对象的元参数 (‘mp’) 项集成。
• 执行MCMC时，应用No-U-Turn（NUTS）核，从SVI首先获得的平均后验估计开始。

估算恒定细胞周期速度的近似点

• 模型凭借求解一阶微分方程$\frac{d}{dt}{s}_g(t)={\beta }_g{u}_g-{\gamma }_g{s}_g$来获得初始洞察，其中${\gamma }_g$是基因依赖的降解率，${\beta }_g$ 和 $u_g$分别是与基因相关的参数。
• 假设未剪接读数$u_g(t)$遵循单谐波周期函数，即$u_g(t)=u_{0g}(1+ϵ\cos (\omega t-{\phi }_{0g}))$，其中 $\omega$表示细胞周期速度，${\phi }_{0g}$是相位偏移，$ϵ$ 是幅度。
• 基于上述假设，剪接读数$s_g(t)$具有相同的函数形式，但幅度和相位经过调整，即$s_g(t)=s_{0g}(1+{ϵ}^{\prime }\cos (\omega t-{\phi }_{ig}))$。其中，${ϵ}^{\prime }=ϵ\cos (\Delta {\phi }_g)$，$\Delta {\phi }_g=({\phi }_g-{\phi }_{0g})$，且 $\tan (\Delta {\phi }_g)=\omega {\gamma }_g^{-1}$。这表明相位差和幅度调整与细胞周期速度$\omega$ 和降解率 ${\gamma }_g$ 相关。
• 假设存在多个条件（或重复实验），且寿命${\tau }_g={\gamma }_g^{-1}$与条件无关，观察到关系${\delta }_{cg}=\tan (\Delta {\phi }_{cg})={\omega }_c{\tau }_g$。这表示相切值${\delta }_{cg}$可看作细胞周期速度${\omega }_c$ 与寿命 ${\tau }_g$ 的乘积。
• 通过奇异值分解 (SVD)，${\delta }_{cg}$许可分解为秩-1 矩阵形式，即${\delta }_{cg}=u_{c}d{v}_g+$更高秩项，其中$u_c$ 和 $v_g$分别是条件和基因的向量，$d$ 是标量。
• 基于 SVD，结果可进一步表达为条件特定的细胞周期速度${\omega }_c$，以逆平均半衰期单位（记为${\omega }_c^{\ast }$）表示，即 ${\omega }_c^{\ast }=u_{c}d{v}_g$。其中 $v_g$是基因的平均值。
• 周期长度以平均半衰期单位表示为$T_c^{\ast }=\frac{2\pi }{{\omega }_c^{\ast }}$，反映了细胞周期的周期性特性。

数据集

• ‘small’：包含 97 个基因。
• ‘medium’：包括 218 个基因。默认使用。
• ‘large’：囊括 1,918 个基因。
• 使用 velocycle.utils.get_cycling_gene_set 函数访问这些人类和老鼠的基因集。