比较基因组学的统计第一类错误

我想梳理两件事:第一,像 RERconverge 这种把“基因速率信号”直接和离散表型做关联的框架,为什么在真实数据里很容易被树上的过程变异误导;第二,以隐藏状态/贝叶斯层级模型为思路的方法,如何把分支背景变异(例如 $\beta_l^2$)放进模型里,从而把这类误导压下去。顺带提醒自己:统计上看到多条分支“同时加速”,并不等于它们因为同一个表型发生了因果趋同。

如果用 RERconverge 这类“速率—性状直接关联”的框架去问:

离散表型(如水生/陆生)是否会在整体上改变某个基因的演化速率?

最容易踩到的坑是 系统发育伪重复(phylogenetic pseudoreplication):祖先分支上只发生过一次的过程变异,被后代整支系继承;到了相关性检验里,却被当成很多次独立证据,再和表型标签对齐,于是 P 值看起来漂亮得不真实。

MuSSCRat(May & Moore 2020)走的是另一条路:把“沿树传递的背景差异”显式写成(祖先或分支上的)$\beta_l^2$(分支背景速率),再和状态效应 $\zeta^2$(离散状态对应的相对速率)一起做后验推断。这样一来:

  • 单系群里“继承下来的一次加速”,更倾向被解释为祖先分支的 $\beta^2$,而不是被写成“表型导致的整体效应($\zeta^2$)”
  • $\zeta^2$ 不再被“同一事件在很多后代里的拷贝”虚假抬高

下面按同一套例子对照:RERconverge 的显著性是怎么在“伪重复 + 过程变异”的组合下被推高的;MuSSCRat 又是怎样靠分支背景项($\beta_l^2$)把背景变动和状态效应分开。


我的思考路径

读 Tribble et al. 2026(An Evolving View of Character Macroevolution)时,我反复看到同一个风险:没有生成模型约束时,把速率信号直接拿去和表型做关联,很容易被 过程变异(process variation) 推出假阳性。我一开始想到用隐藏状态去吸收未知变动,但 RERconverge 这种回归式流程很难自然加入隐藏状态。

为了找一个更具体的对照例子,我转去看 May & Moore 2020 的 MuSSCRat。它用贝叶斯层级模型评估离散特征对连续速率的影响,同时给每条分支一个背景项来接住背景变动。真正理解机制时我曾被 CTMC 这类术语卡住,后来换了入口:不追推导,从数据生成去理解参数——总方差可以拆成状态效应 $\zeta^2$ 和分支背景 $\beta^2$。把 RERconvergeMuSSCRat 放在同一数据结构里对照,假阳性更像是在这组条件下出现:单系群的过程变异没有被隔离,而检验阶段又弱化了树的非独立性;MuSSCRat 的联合拟合则能把“祖先分支的一次 $\beta^2$ 跃升”从“水生的整体效应”里拆出来。


直接关联 vs 隐藏状态/层级建模

数据结构

  • 离散特征 X:栖息地(0 = 陆生,1 = 水生)
  • 连续响应 Y:某个基因在各分支上的相对演化速率(RER 或等价的“分支速率信号”)

例子

鲸豚类(海豚、鲸鱼等多个物种)在 OR9Q2 上普遍呈现高 RER;海牛、海豹等其他水生谱系并不呈现同样的加速模式。

两种方法的路径差异

  • **RERconverge**:先用全基因组平均枝长做代数标准化,再把标准化后的 Y 与 X 做相关性检验。
  • **MuSSCRat**:不把背景变动交给预处理,而是在概率模型里为每条分支引入 $\beta_l^2$,并与 $\zeta^2$ 联合估计。

假阳性是怎么被制造出来的:系统发育伪重复与单系群过程变异

当单系群(clade)内部有一段没被隔离的过程变异时,RERconverge 这类“速率—性状直接关联”的框架就很容易在真实数据里给出夸张的显著性。

  • 真实历史里发生了什么:假设在鲸豚类共同祖先分支上,因为一次全基因组事件或强烈瓶颈,OR9Q2 出现了一次长期的速率跃升,然后被整个鲸豚支系继承。它在演化史上只发生过一次,也可能和“进入水生环境”没有因果关系。
  • **RERconverge 的统计视角:相关性检验并不认识“继承事件”。它把多个鲸豚物种的叶节点当作相互独立的点,于是会看到很多 **(X=水生, Y=RER 高) 的配对,秩相关自然显著。
  • 为什么这会把因果说歪:要说“水生驱动加速”,你更希望看到多次独立重复,比如几个互不相近的水生类群都各自出现加速。若现实里主要是鲸豚高、海牛海豹不高,更贴近数据的说法往往只是“OR9Q2 在鲸豚里加速了”,而不是“水生普遍带来加速”。
  • 盲点出在哪里:当速率信号被标准化/汇总成可回归的数值后,树的拓扑与非独立性在检验阶段很容易被弱化,从而难以区分:
    • 一次祖先事件在多个后代中的继承(伪重复),与
    • 多次独立演化中的趋同(更接近独立证据)。

麻烦的不是“随机噪声”,而是演化历史的结构没有写进检验假设:过程变异和伪重复叠在一起,显著性就会被推高。


MuSSCRat:把“背景方差”从统计检验里拿出来,交给模型参数

面对“鲸豚类整支系一致加速、但海豹/海牛并不如此”的数据,MuSSCRat 不需要把它写成“水生效应”。它把连续性状(或连续速率信号)的演化写成“方差怎么生成”的问题,并把总方差拆成两类因子:

\[\text{TotalVariance}_l = \beta_l^2 \times \zeta_{X(l)}^2\]

其中:

  • $\zeta^2$:离散状态对应的相对速率(状态 0 vs 状态 1 的差异)
  • $\beta_l^2$:分支背景速率(每条分支自己的背景变动,相当于过程变异的“容器”)

MuSSCRat 里,“鲸豚类后代普遍高”不必被解释成“水生整体更快($\zeta_1^2$ 很大)”。模型允许你说:

水生整体未必更快;高方差来自鲸豚类祖先分支的一次局部背景跃升,并沿树向下游传递。

为什么模型更愿意抬“祖先分支的 $\beta^2$”

在后验拟合里,模型发现鲸豚类内部许多叶节点的演化速率信号都偏高时,典型权衡是:

  • 路线 A:抬高 $\zeta_1^2$(解释成“水生整体更快”)
  • 路线 B:只抬高鲸豚类共同祖先那条分支的 $\beta^2$(解释成“该祖先分支上的一次背景跃升”,并通过沿树传递解释整支鲸豚)

如果走路线 A,会连带要求所有水生谱系都更快,于是海豹、海牛这些同样标注为水生、但速率不高的分支会出现大面积错配(惩罚很大)。

相反,走路线 B 更省:只调整一个祖先分支的 $\beta^2$,就能解释整条鲸豚类支系的高信号,同时不需要海豹/海牛去拟合同样的加速。

后验往往更倾向路线 B:把“单系群继承”归因到祖先分支的 $\beta^2$,让全局的 $\zeta^2$ 不至于被伪重复堆高。

输出怎么读更稳妥

与其把“秩相关/秩和检验的 P 值”当作结论,不如直接看状态效应参数是否被数据支持。P 值很容易被单系群内整齐一致的信号推到极显著,但那未必对应你想要的解释。更稳妥的读法是:

  • $\zeta_1^2$ 与 $\zeta_0^2$ 的后验分布是否明显分离
  • 或者用贝叶斯因子比较“有状态效应”与“无状态效应”的模型支持度

术语与必要的数学思维(卡壳时再读)

这一节不做推导,只把术语翻译成“它在算什么”。主线不依赖这里。

1) CTMC:不是生物学理论,而是“计算概率的引擎”

文献里说一个模型在 CTMC(连续时间马尔可夫链)框架下,通常只意味着:

给定“状态”与“瞬时转移速率”,算法能在任意枝长 $t$ 上计算从起点状态到终点状态的概率。

CTMC 需要三个前提(记住它们就够用):

  1. 离散状态:对象是有限分类(例如 A/T/C/G,或表型 0/1)。
  2. 连续时间:系统发育树的枝长是时间,状态改变可以发生在任意瞬间。
  3. 无后效性(马尔可夫性):未来只依赖“当前状态 + 固定速率”,不依赖更久远的历史细节。

最常见的写法是一个速率矩阵 $Q$。以二元表型 $0,1$ 为例:

$$ Q = \begin{bmatrix}

  • & q_{01} q_{10} & - \end{bmatrix} $$
  • $q_{01}$:从 0 到 1 的瞬时转移速率
  • $q_{10}$:从 1 到 0 的瞬时转移速率
  • 对角线:定义为该行非对角元素之和的负数(保证概率守恒)

理解到这里就够了:CTMC 负责把 $Q$ 和枝长 $t$ 变成“这条分支上状态如何变化”的概率;更深的矩阵指数/微积分细节,不是理解本文主线所必需的。

2) 边缘化(Marginalization):把不可观测的维度“加总消掉”

隐藏状态模型里只能观测到 $O$,但模型内部会引入不可直接观测的 $H$。

如果底层计算的是联合概率 $P(O,H)$,你想要观测层面的概率 $P(O)$,就要把所有可能的 $H$ 加总:

\[P(O) = \sum_H P(O,H)\]

直觉翻译:我不知道隐藏状态是哪一个,但我可以把“所有可能历史”的概率都算进去,再加起来。

3) 直接关联 vs 层级建模:差别不在“更复杂”,而在“背景异质性有没有显式进模型”

“直接做关联”的思路可以被理解成:先把树与速率信号处理成一个可用于检验的观测量,再在其上做条件推断。

例如(只看结构):

\[P(M \mid O, \hat{T}) \propto P(O \mid M, \hat{T})P(M)\]

层级/隐藏状态类模型则是在概率模型内部显式加入“看不见但客观存在”的异质性维度(例如分支背景项 $\beta_l^2$),并在后验里让数据去分配解释权:哪些信号应该归到“状态效应”,哪些信号更像“分支背景异质性”。

4) MCMC:不是“求一个最优解”,而是“在后验里采样”

当参数很多、解析解很难写出来时,MCMC 用随机游走的方式在后验分布里采样。最终得到的不是一个点估计,而是一整个分布(这也是为什么 MuSSCRat 更自然地表达不确定性)。


收束:频率 vs 贝叶斯(以及这篇笔记要保留的警报器)

两种范式在根本上回答不同的问题:

  1. 频率学派常见形式:$P(\text{Data}\mid H_0)$ —— 解释的是:如果零假设为真,看到当前数据(或更极端)的概率多大。
  2. 贝叶斯常见形式:$P(H\mid \text{Data})$ —— 解释的是:看到当前数据后,你的假设为真的概率分布是什么。

问题通常不在工具,而在误用:把 P 值当作“假说为真的概率”是一种常见的逻辑跳跃。

做“离散表型 vs 连续速率”的比较分析时,如果看到以下症状,就要警惕“速率—性状直接关联”得到的显著性可能是伪重复撑起来的:

  • 显著性由某个单系群内部“整齐一致”的信号驱动(多个物种同高/同低,像同一祖先事件的后遗症),而不是多个彼此独立的水生类群各自独立复现同一模式
  • 把离散状态当作独立重复:表型标签在多个物种上相同,但速率信号在系统发育上高度相关(伪重复),相关性检验却把它们当独立点

这往往说明你面对的是 过程变异 + 树上的非独立性:更合适的处理方式是把它显式写进模型(例如用祖先/分支上的 $\beta_l^2$ 吸收单系群背景),而不是交给“标准化 + 丢开拓扑的相关性检验”硬解释。


另一类“看似假阳性”的强信号:统计学并发加速 ≠ 生物学因果趋同

Tribble et al. 2026 提到“缺少联合建模会放大过程变异”的时候,我还想到另一种更常见的错位:贝叶斯竞争模型(例如 PhyloAcc 这类把前景分支组合当成候选并比较 $H_0/H_1/H_2$ 的框架)在某些数据格局下会给出很强的统计支持,但统计支持不等于因果解释成立。

统计学输出与生物学因果的错位

如果随意挑选的一组前景(例如“奶牛 + 蝙蝠”)在真实数据里确实同时发生了显著的速率加速,那么模型输出高后验概率与极高贝叶斯因子去支持“多前景独立加速”的 $H_2$,在统计学上完全合理:模型回答的就是“这些分支是否同时出现了独立的速率改变?”而不是“它们是否共享某个离散表型的同一因果机制?”。

错位发生在解释阶段:把“统计学上并发的速率改变(parallel evolution in rates)”直接说成“共享离散表型导致的趋同适应”。当某个前景根本不具备你假设的表型(例如把“回声定位”套到奶牛身上)时,因果链条已经断了,但模型只看序列,不可能替你识别这种语义冲突。

$\beta_l^2$ 这类分支背景参数能压住“只有单条分支加速、却被写成表型效应”的问题,因为信号可以被吸收到某条分支的背景项里;但它挡不住另一种情况:多条分支确实同时加速,只是动因不一样。对只看序列的速率模型来说,那仍然是“真实存在的并行信号”。

无生物学关联的速率趋同常见吗

在全基因组尺度上,两条表型差异很大的谱系在同一个基因上出现并行加速并不罕见;很多时候,它来自与“宏观表型相似”无关、但会产生强速率信号的机制。

  • 净化选择的平行放松(parallel relaxation of purifying selection)与假基因化:当基因功能变得冗余,维持保守性的净化选择压力减弱,突变迅速积累,表现为速率升高。这个机制与“水生/嗅觉受体”语境天然兼容:水生环境下嗅觉或某些受体通路功能边际下降时,相关基因更容易出现退化式加速。此时“多个水生类群都加速”也未必是同一正向选择驱动的趋同适应,而可能是“功能不再必要”导致的并行退化。
  • 基因多效性带来的不同选择压力:同一基因参与多个生理过程时,不同谱系可能因为完全不同的生态压力在不同结构域上承受选择,但在“全长速率”统计量里都被记录为“加速”。把这种并发加速直接归为同一离散表型的因果效应,风险很高。
  • 基因组结构层面的中性偏置(例如 gBGC):在重组率高的区域,GC 偏向性基因转换会推动固定突变累积,使某些区域在速率模型里呈现出并行加速;这一类信号可以很强,但与关注的表型无关联。

面对真实存在的随机并行加速,应该怎么处理

检测到“强支持 $H_2$”之后,更稳妥的做法是把“统计学并发”和“生物学因果”拆开验证。

1) 经验零分布(permulations 思路):如果某个前景组合(例如蝙蝠 + 海豚)在 PhyloAcc 输出极高 BF,一个自然的问题是:这种程度的并发加速,是否显著超出“该基因在全树上随机出现并发加速”的本底概率?做法是:在物种树上随机生成大量与真实前景在拓扑性质上匹配的伪前景组合(例如保持前景数量与拓扑距离),对同一基因重复计算 BF,得到经验零分布。只有当真实前景的 BF 明显高于该零分布(例如超过 95% 分位数),才更有把握拒绝“随机并行加速”这一类零解释。

2) 拆解速率来源:正向选择 vs 功能丧失:仅知道“RER 升高/速率加快”并不足以支持适应性趋同。一个常用的后验鉴别是用密码子模型(例如 PAML 的分支-位点模型)去看加速来自哪里:$\omega$ 若趋近 1 并伴随提前终止密码子,更像平行假基因化(功能丧失);$\omega>1$ 才更像正向选择。即便 $\omega>1$,也需要进一步看加速位点是否集中在相同功能结构域:如果不同前景的加速位点落在不同结构域,更像多效性下的独立事件,而不是同一离散表型驱动的同一因果路径。