稳态模型：RNA velocity的数学基础
好，我们正式进入RNA velocity的数学核心。说实话，很多初学者看到微分方程就头疼。但别担心，这一章我会带着你把这个模型拆开揉碎，讲明白它到底在算什么。
RNA velocity之所以能预测细胞状态的变化方向，核心就是靠一个描述基因转录与降解动力学的常微分方程组。我最早看到这个模型时也觉得有点抽象，后来自己用模拟数据跑了一遍，才真正理解它的美妙之处。
3.1 核心微分方程
我们先从最基础的模型说起。对于一个给定的基因，它有两个状态变量：
- u(t)：未剪接（unspliced）mRNA的数量，代表刚刚转录出来的前体RNA
- s(t)：已剪接（spliced）mRNA的数量，代表成熟的功能性RNA
它们的变化速率由以下方程组描述：
du/dt = α(t) - β · u(t)
ds/dt = β · u(t) - γ · s(t)这里三个参数是关键：
| 参数 | 含义 | 我的经验值 |
|---|---|---|
| α(t) | 转录速率（可能随时间变化） | 不同基因差异巨大，范围可跨越3个数量级 |
| β | 剪接速率常数 | 通常假设为常数，我个人常设1.0 |
| γ | 降解速率常数 | 这是需要估计的核心参数 |
为什么是这两个方程？很简单：
- 转录产生未剪接RNA，所以du/dt的第一项是+α(t)
- 未剪接RNA通过剪接变成成熟RNA，这同时减少u、增加s
- 成熟RNA会被降解，所以s项有一个-γ·s的损失
我个人习惯把这两个方程想象成一个“水池模型”：u是上游水库，s是下游水库。α是上游的来水，β是闸门开度，γ是下游的渗漏。你想想看，只要监测两个水库的水位变化，就能推断出闸门和渗漏的情况——RNA velocity也是同样的逻辑。
3.2 稳态假设：让问题可解
好，方程有了。但问题来了：我们只有一次snRNA-seq的快照数据，每个细胞只能测到u和s的当前值，根本不知道它们随时间的变化轨迹。怎么办？
这时候就需要引入稳态假设。这个假设说的是：在足够长的时间尺度上，基因的表达水平会达到一个稳态。在这个稳态下，转录、剪接、降解达到平衡，u和s都不再变化。
数学上，稳态意味着：
du/dt = 0, ds/dt = 0代入原方程得到：
0 = α - β · u
0 = β · u - γ · s解这个方程组，我们会得到一个关键的关系：
u = α / β
s = (β / γ) · u说白了，在稳态下，s与u成正比。这个比例系数就是β/γ。
我在项目中遇到过这样的情况：有些基因在特定细胞类型中永远达不到稳态，比如细胞周期基因。这种情况下强行用稳态模型就会出问题。所以后来有了动态模型（我们后面会讲），但稳态模型仍然是理解和入门的基础。
3.3 从数据中估计速率
现在我们知道，稳态下s和u应该是一条过原点的直线。但实际数据中，每个细胞的s和u值会散落在直线周围。那些偏离直线的点，就是“正在发生变化”的细胞。
具体怎么算？看这个流程：
对于每个基因：
1. 收集所有细胞的u和s值
2. 取u值最高的前X%的细胞（通常是前5%）做线性回归
3. 回归线的斜率 = β/γ
4. 设定β = 1（或可估计），计算γ = β / 斜率
5. 对每个细胞计算残差：δ = s - (β/γ) · u
6. velocity = β · u - γ · s = -γ · δ这里的核心逻辑是：高表达量的细胞更可能处于稳态，因为它们有更多的机会达到平衡。所以我一般取前5%的高表达细胞来拟合回归线——这是La Manno等人论文中的标准做法。
关键结论：当δ > 0时，s偏高，说明剪接已经发生但降解还没跟上，velocity为负（正在减少）。当δ < 0时，s偏低，说明新转录的RNA还在路上，velocity为正（正在增加）。
3.4 稳态模型的局限与我的避坑指南
这里我必须泼一盆冷水：稳态模型虽然优雅，但并不完美。我从业这些年，踩过的坑主要包括：
- 假设所有基因都能达到稳态——实际上很多基因的转录是脉冲式的，比如应激反应基因
- 所有细胞共享同一套β和γ——这个假设太强了
- 对低表达基因不友好——表达量低的基因信噪比很差
我曾经踩过一个坑：拿一个在T细胞中表达很低的趋化因子基因去跑velocity，结果出来的箭头全是乱的。后来发现就是这个基因的表达量太低了，稳态回归的斜率根本不可靠。从那以后，我给自己定了个规矩：对于平均表达量低于0.1的基因，直接过滤掉。
避坑指南：不要直接拿原始counts做线性回归。我个人习惯先做log1p变换，然后把表达量按基因归一化到[0,1]区间。这样能减少不同基因间表达量级差异带来的偏差。
另外还有个重要细节：我们前面说设定β=1。这其实只是给定了时间尺度。β的具体值会影响velocity的绝对值，但不影响方向。而RNA velocity最核心的用途是判断方向，所以这通常不是问题。
我的技巧：在实际项目中，我经常做一个小验证——把所有基因的velocity加起来画一个总体的流场图，看看是不是和已知的发育轨迹一致。如果发现某个已知的终末分化细胞群却指向未分化状态，那就说明模型可能出了问题，需要检查那个通路的基因表达是否正常。
3.5 稳态模型知识框架
好的，让我们用一张SVG图把整个稳态模型的知识体系串起来。这张图展示了从原始测序数据到最终RNA velocity计算的核心逻辑。
这张图把整个流程串起来了。你可以看到，从原始数据出发，经过预处理和稳态假设，再到线性回归和残差分析，最后计算出每个细胞每个基因的velocity。这个流程看起来简单，但每个环节都有不少需要注意的细节。
我自己最初学习时，是先把这张图手画了一遍，记住每个箭头表示的含义。然后对照着实际的代码一行一行看，把数学和工程对应起来。嗯，这个方法你可以试试，确实能帮助理解。
最后说一句——稳态模型是整个RNA velocity领域的基石。后面所有的改进模型（动态模型、EM模型等）都是在它的基础上做扩展。你把这章吃透了，后面30章都会轻松许多。