【论文阅读笔记】VeloCycle

流形

• 通过单细胞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({(Sc,Uc)}∣θ)P(\left\{(S_c, U_c)\right\} | \theta)P({(Sc​,Uc​)}∣θ) 由以下表达式组成：
  P({(Sc,Uc)}∣θ)=∏gcP(Sgc,Ugc∣ω(φ),φc,vg,βg,γg,αg)P(Sgc,Ugc∣θ)=Ps(Sgc∣vg,αgs,φc)×Pu(Ugc∣ω(φ),βg,γg,vg,φc,αgu)Ps(Sgc∣…)=NB(sg(φc)=F[vg,φc],αgs)P(Ugc∣…)=NB(ug(φc)=G[ω(φc),βg,γg,vg,φc],αgu)P(\left\{(S_c, U_c)\right\} | \theta) = \prod_{gc} P(S_{gc}, U_{gc} | \omega(\varphi), \varphi_c, v_g, \beta_g, \gamma_g, \alpha_g) \\ P(S_{gc}, U_{gc} | \theta) = P_s(S_{gc} | v_g, \alpha_g^s, \varphi_c) \times P_u(U_{gc} | \omega(\varphi), \beta_g, \gamma_g, v_g, \varphi_c, \alpha_g^u) \\ P_s(S_{gc} | \ldots) = \text{NB}(s_g(\varphi_c)=F[v_g, \varphi_c], \alpha_g^s) \\ P(U_{gc} | \ldots) = \text{NB}(u_g(\varphi_c)=G[\omega(\varphi_c), \beta_g, \gamma_g, v_g, \varphi_c], \alpha_g^u) P({(Sc​,Uc​)}∣θ)=gc∏​P(Sgc​,Ugc​∣ω(φ),φc​,vg​,βg​,γg​,αg​)P(Sgc​,Ugc​∣θ)=Ps​(Sgc​∣vg​,αgs​,φc​)×Pu​(Ugc​∣ω(φ),βg​,γg​,vg​,φc​,αgu​)Ps​(Sgc​∣…)=NB(sg​(φc​)=F[vg​,φc​],αgs​)P(Ugc​∣…)=NB(ug​(φc​)=G[ω(φc​),βg​,γg​,vg​,φc​],αgu​)
  其中$\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ωt,{φc},{vgt},{βg},{γg},{αg})P(\{v\omega_{ t}, \{\varphi_c\}, \{v_{gt}\}, \{\beta_g\}, \{\gamma_g\}, \{\alpha_g\})P({vωt​,{φc​},{vgt​},{βg​},{γg​},{αg​}) 被分解为多个独立分量的乘积，其形式为：
  P({vωt,{φc},{vgt},{βg},{γg},{αg})=∏cP(vωt)P(φc)P(vgt)P(βg)P(γg)P(αg)P(\{v\omega_{ t}, \{\varphi_c\}, \{v_{gt}\}, \{\beta_g\}, \{\gamma_g\}, \{\alpha_g\}) = \prod_c P(v\omega_{ t}) P(\varphi_c)P(v_{gt}) P(\beta_g)P(\gamma_g) P(\alpha_g) P({vωt​,{φc​},{vgt​},{βg​},{γg​},{αg​})=c∏​P(vωt​)P(φc​)P(vgt​)P(βg​)P(γg​)P(α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$ 之间的关系。
• 后验因子分解如下：
  P({νωt},{φc},{νgt},{βg},{γg},{αg})=P({γg},{νωt})∏gP(βg∣γg)P(αg)∏tP(νωt)P(νgt)∏cP(φc)\begin{align*} P\left(\{\nu \omega_t\}, \{\varphi_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(\varphi_c) \end{align*} P({νωt​},{φc​},{νgt​},{βg​},{γg​},{αg​})​=P({γg​},{νωt​})g∏​P(βg​∣γg​)P(αg​)t∏​P(νωt​)P(νgt​)c∏​P(φc​)​
• 具体公式是：
  x≡[log⁡(γ1),log⁡(γ2),…,log⁡(γng),νω0,νω1,…,νωnt−1]Σ=F^F^⊤+diag(d^)where F^∈R(ng+nt)×k,with k=5P({log⁡(γg)},{νωt})=P(x)=MultivariateNormal(m^,Σ)μlog⁡βg∣γ=μ^log⁡βg+ρ^g⋅μ^log⁡βg⋅log⁡(γg)−μ^log⁡γgσlog⁡γgwith ρ^g∈[0,1]σlog⁡βg∣γ=μ^log⁡βg1−ρ^g2P(log⁡(βg)∣log⁡(γg))=N(μlog⁡βg∣γ,σlog⁡βg∣γ2)P(φc)=N([φ^xc,φ^yc],[1,1])P(νgt)=N(μ^gtν,σ^gtν2)P(αg)=Delta(αg^)\mathbf{x} \equiv \left[ \log(\gamma_1), \log(\gamma_2), \ldots, \log(\gamma_{n_g}), \nu \omega_0, \nu \omega_1, \ldots, \nu \omega_{n_t-1} \right] \\ \boldsymbol{\Sigma} = \hat{\mathbf{F}} \hat{\mathbf{F}}^\top + \mathrm{diag}(\hat{\mathbf{d}}) \quad \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}}, \boldsymbol{\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} } \quad \text{with } \hat{\rho}_g \in [0,1] \\ \sigma_{\log \beta_g | \gamma} = \widehat{\mu}_{\log \beta_g} \sqrt{1 - \widehat{\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(\varphi_c) = \mathcal{N}([\widehat{\varphi} \mathbf{x}_c, \widehat{\varphi} y_c], [1,1]) \\ P(\nu_{gt}) = \mathcal{N}(\widehat{\mu}_{gt}^{\nu}, \widehat{\sigma}_{gt}^{\nu 2}) \\ P(\alpha_g) = \text{Delta}(\widehat{\alpha_g}) x≡[log(γ1​),log(γ2​),…,log(γng​​),νω0​,νω1​,…,νωnt​−1​]Σ=F^F^⊤+diag(d^)where F^∈R(ng​+nt​)×k, with k=5P({log(γg​)},{νωt​})=P(x)=MultivariateNormal(m^,Σ)μlogβg​∣γ​=μ^​logβg​​+ρ^​g​⋅μ^​logβg​​⋅σlogγg​​log(γg​)−μ^​logγg​​​with ρ^​g​∈[0,1]σlogβg​∣γ​=μ​logβg​​1−ρ​g2​​P(log(βg​)∣log(γg​))=N(μlogβg​∣γ​,σlogβg​∣γ2​)P(φc​)=N([φ​xc​,φ​yc​],[1,1])P(νgt​)=N(μ​gtν​,σgtν2​)P(αg​)=Delta(αg​​)

模型实现

• 模型实现旨在估算联合后验概率分布的近似值，涉及角细胞周期速度 ($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 函数访问这些人类和老鼠的基因集。