SNeP估算Ne历史曲线的原理

Posted on | Post modified | In

简单来说，这个软件的原理是：两个位点之间的“连锁不平衡”程度，就像是记录它们在过去某段时间里共同历史的“快照”。通过分析不同距离的位点对的LD程度，我们就能像拼拼图一样，重建出不同历史时期的有效群体大小。

为什么LD能反映Ne？

核心概念辨析：LD的本质是什么？—— 不是物理连锁，而是统计关联

  • 物理连锁 是一个固定的、结构性的概念。它指两个位点在同一条染色体上，且物理位置靠得很近。只要它们没被重组事件分开，它们就始终“连锁”在一起。这是染色体的一种物理状态。
  • 连锁不平衡（LD） 是一个动态的、统计性的概念。它衡量的是两个位点的等位基因在种群中共同出现的非随机性。它描述的是基因频率分布的模式。简单说，就是两个在染色体上位置相近的遗传位点（比如两个SNP），它们的等位基因不是随机组合的。例如，如果位点A的“A”等位基因总是和位点B的“T”等位基因一起出现，而很少和“C”一起出现，那么这两个位点就处于LD状态。 关键区别：即使两个位点紧密连锁（物理上很近），如果它们的等位基因是随机组合的（比如A总是50%概率和T在一起，50%和C在一起），那么它们的LD值也为0，即处于“连锁平衡”状态。

LD如何产生和衰减？

  • 产生：一次新的突变产生时，它必然与它所在那条染色体上周围的所有位点处于完全LD状态。
  • 衰减：主要靠两个过程：
    • 重组：在减数分裂中，染色体之间的交叉互换（重组）会打断位点之间的联系。两个位点物理距离越远，发生重组把它们分开的概率就越大。
    • 遗传漂变：在小群体中，等位基因频率的随机波动也会改变它们之间的关联模式。群体越小（Ne越小），遗传漂变的力量越强。

关键公式（Hill, 1981）与“LD期望值（E(r²)）”

理论推导出一个简洁的关系：在一个理想的随机交配群体中，两位点间的LD期望值（通常用 r² 衡量）与过去的有效群体大小（Ne） 和重组率（c） 相关： E(r²) ≈ 1 / (1 + 4Ne * c) 这里 c 就是两位点间的重组率，与物理距离成正比。

如何理解 “LD期望值（E(r²)）”？ 需要澄清的是，LD期望值（E(r²)）不等于两个位点“物理联锁”的概率。 E(r²) 中的 E 代表 期望，这是一个概率论和统计学概念。它指的是：在一个符合特定理论假设（如随机交配、无选择、恒定群体大小等）的理想群体中，对于给定重组率 c 的两个位点，我们预期能观察到的平均关联强度（r²值）是多少。 可以这样类比：

  • 抛一次硬币：得到正面或反面，这是一次观测结果。
  • 抛无穷多次硬币：正面向上的比例期望值是50%。这个“50%”就是期望值。
  • 同理，在真实数据中，我们观测到的是成千上万对特定距离的SNP的具体r²值（这些值有高有低）。而 E(r²) 是理论告诉我们，在理想情况下，这些r²值应该围绕哪个平均水平上下波动。 所以，E(r²) 是关于“等位基因关联强度”的理论预期值，而不是“物理联锁的概率”。

为什么 E(r²) 与 Ne 和 c 有关？—— 一个动力学比喻 想象一个住满人的（Ne）大广场，每个人手上都拿着两根长度不同的绳子（代表染色体上两个位点）：

  • 短绳子（c小）： 两根绳子离得很近，很容易被绑在一起或保持在一起。
  • 长绳子（c大）： 两根绳子离得远，很容易被分开。
  • 广场上的人（Ne）： 人越多，环境越复杂，随机走动时绳子被意外缠住或解开的方式就越多。 E(r²) 衡量的是，在广场上随机找两个人，观察他们手上那两根绳子被“绑在一起”的比例的期望值。
  • 如果广场人很少（Ne小）：随机性（漂变）很强。可能纯粹因为运气，大部分人的短绳子都刚好被同一种方式绑住了，那么你观察到的“绑在一起”的比例（r²）就会很高。E(r²) 理论值也就高。
  • 如果广场人很多（Ne大）：随机性很弱。绳子是否绑在一起主要看它们本身的特性（距离c）。长绳子很难绑在一起，所以观察到的比例（r²）就很低。E(r²) 理论值也就低。
  • 重组率c： 它直接决定了绳子“被解开”的难易程度。c越大（绳子越长/距离越远），天生就越容易被解开，所以E(r²)的期望值就越低。

通过以上的动力学过程，公式 E(r²) ≈ 1 / (1 + 4Ne * c) 进行了数学上的精炼表达。它的直觉是：如果历史上Ne很大，遗传漂变弱，打破LD主要靠重组。那么，对于同样距离（c值相同）的位点对，它们在过去经历的重组次数相对固定，所以LD水平（r²）会较低。反之，如果历史上Ne很小，遗传漂变很强，它会更快地改变等位基因组合，使得LD水平（r²）较高。

SNeP如何操作 —— 从数据到历史

软件的核心算法就是基于上述理论，但需要处理真实数据的复杂性。过程如下：

  1. 按距离分组：
    • 软件计算全基因组上数百万对SNP之间的物理距离和LD值（r²）。
    • 它不会一对一对地看，而是将所有SNP对按物理距离（比如0-0.1 cM, 0.1-0.2 cM…）分组。每个距离区间代表一个特定历史时期。
      • 为什么？ 因为距离≈重组率≈时间。两个位点越近，在过去发生重组从而断开联系所需的世代数就越多（即它们反映的是更近的历史）。位点越远，重组更容易发生，它们反映的是更久远的历史。
  2. 为每个距离区间计算一个“平均r²”：
    • 对于每个距离区间（比如代表过去~50世代的区间），软件计算该区间内所有SNP对的平均r²值。这个“平均r²”就被视为该历史时期LD水平的一个估计。
  3. 代入公式，反向求解Ne：
    • 对于第 i 个距离区间，软件知道：
      • r²_i（观察到的平均LD）
      • c_i（该区间中心点对应的重组率，由物理距离通过映射函数换算而来）
    • 将 r²_i 和 c_i 代入修改后的理论公式（经过多种校正，见下文），就可以反解出 Ne(t)。这里的 t 就是这个距离区间所代表的历史世代数，通常用 t = 1/(2c) 来近似估计（意思是大约从 t 代以前开始，重组有显著机会打断这两个位点）。
  4. 绘制Ne随时间变化的曲线：
    • 对每一个距离区间重复上述步骤，我们就得到了一系列点：(世代数t, 有效群体大小Ne(t))。
    • 将这些点连接起来，就形成了一条 “Ne历史变化曲线”，X轴是回溯的世代数，Y轴是该时期的估计有效群体大小。

重要的校正（SNeP的改进之处）

原始理论公式很简洁，但现实数据有偏差，所以SNeP引入了关键校正：

  1. 样本量校正：用小样本计算的r²会系统性偏高。软件使用 r²_adjusted = r²_observed - 1/(2 * 样本数) 进行校正。
  2. 突变校正：古老的突变会降低观察到的LD。软件引入了一个参数α来校正。
  3. 单倍型相位不确定性校正：当基因型数据未被定相（即不知道等位基因在同一条染色体上的排列）时，LD估计会偏低。软件对此进行补偿。
  4. 灵活的分箱策略：确保每个时间区间内有足够多的SNP对，使平均r²更稳定。

总结与直觉比喻

整个过程好比是：

  • 基因组是一本厚重的、由无数个字母（SNP）写成的历史书。
  • LD（r²） 是测量相邻两个字母在抄写传承过程中，保持原样的“关联强度”。
  • 距离决定了你读的是书的哪一章：读紧密相邻的字母（距离近），你看到的是最近几页（近期历史）的抄写错误和关联；读隔得很远的字母（距离远），你看到的是很久以前章节（远古历史）的关联。
  • SNeP 就像一个聪明的统计学家，它测量整本书中所有不同距离的字母对之间的平均关联强度，然后利用一个已知的物理定律（重组和漂变如何影响关联），反推出在书写这本书的不同历史时期，抄写员团队的大小（Ne）是如何变化的。

对于SNeP软件而言：

  1. 它观测：全基因组上数百万对SNP在当前群体中的实际LD水平（观测到的r²）。
  2. 它利用理论：这个理论（公式）描述了在过去特定历史条件下（由Ne和c决定），LD的期望水平 E(r²) 应该是多少。
  3. 它反向求解：既然我现在知道了结果（观测r²），也知道了距离c，我就能倒推出导致这个结果的历史原因之一——过去的有效群体大小Ne。 简而言之：LD期望值 E(r²) 不是物理概率，而是一个将“群体历史参数（Ne）”与“可观测的遗传统计量（r²）”连接起来的理论桥梁。SNeP正是走过了这座桥，从我们今日所见的数据，回溯了种群的历史。

最重要的一点是：这个方法估算的远古Ne（曲线的左端）和近期Ne（曲线的右端）误差较大，但整个曲线的变化趋势（是上升、下降还是瓶颈）通常是可靠且有极高生物学意义的。这正是该工具的强大之处。

Ref：https://www.frontiersin.org/journals/genetics/articles/10.3389/fgene.2015.00109/full