跳转到内容

连锁不平衡

快速概览

Hardy–Weinberg 平衡描述了单个位点的稳定性,而连锁不平衡(Linkage Disequilibrium, LD)则描述了不同位点之间的非随机关联。LD 是理解基因组结构、单倍型块以及 GWAS 信号来源的核心概念。

  • 掌握 LD 的数学度量:$D$, $D'$ 和 $r^2$
  • 理解重组(Recombination) 如何作为"打散"LD 的核心力量
  • 掌握单倍型块(Haplotype Blocks) 的块状景观直觉
  • 了解 LD 在基因型插补(Imputation) 和精细定位中的关键作用
所属板块 分析方向与案例

把基础对象与算法方法重新放回真实分析任务与工作流。

阅读目标 帮助建立阅读上下文

先判断这页与你当前问题的关系,再决定是否深入展开。

建议前置 先建立相关基础对象与方法直觉

建议先建立相关基础对象与方法直觉,再进入本页。

如果两个位点是独立遗传的,那么它们同时出现某种等位基因组合的频率应该等于各自频率的乘积。

连锁不平衡(LD) 是指在群体中,不同位点的等位基因之间存在非随机的统计相关性

  • 直觉:如果你知道位点 A 是 G,你就能以很高的概率猜到位点 B 是 A。这意味着这两个位点”连在一起”了。

考虑两个双等位基因位点:位点 1 的等位基因为 A/aA/a,频率为 pAp_Apap_a;位点 2 的等位基因为 B/bB/b,频率为 pBp_Bpbp_b。两个位点共有四种可能的单倍型(Haplotypes)

单倍型观察频率期望频率(独立假设下)
ABABPABP_{AB}pApBp_A \cdot p_B
AbAbPAbP_{Ab}pApbp_A \cdot p_b
aBaBPaBP_{aB}papBp_a \cdot p_B
ababPabP_{ab}papbp_a \cdot p_b

如果两个位点完全独立(连锁平衡),则 PAB=pApBP_{AB} = p_A p_B。任何偏离都意味着存在连锁不平衡。

LD 系数 $D$
观察到的单倍型频率与期望频率之差。$D = P_{AB} - p_A p_B$。
标准化系数 $D'$
将 $D$ 按其理论最大值缩放。$|D'| = 1$ 表示完全连锁(即四个可能的单倍型中至少有一个从未出现过)。
相关系数 $r^2$
GWAS 中最常用的指标。反映了一个位点对另一个位点的预测能力。$r^2 = 1$ 表示两位点信息完全冗余。

DD 是最基本的 LD 度量,但它有几个重要的局限性:

D=PABpApB=PAbpApb=(PaBpapB)=(Pabpapb)D = P_{AB} - p_A p_B = P_{Ab} - p_A p_b = -(P_{aB} - p_a p_B) = -(P_{ab} - p_a p_b)

  • 取值范围不固定DD 的取值范围取决于等位基因频率。当 pAp_ApBp_B 接近 0 时,DD 的最大可能值趋近于 0,这导致不同位点之间的 DD 值不可直接比较。
  • 对称性DD 的符号取决于等位基因的命名方式,互换 A/aA/a 会改变 DD 的符号。

为了解决 DD 的取值范围问题,引入标准化系数:

D=DDmaxD' = \frac{D}{D_{\max}}

其中 DmaxD_{\max} 是在给定等位基因频率下 DD 的理论最大值:

Dmax={min(pApb,papB)if D>0min(pApB,papb)if D<0D_{\max} = \begin{cases} \min(p_A p_b, p_a p_B) & \text{if } D > 0 \\ \min(p_A p_B, p_a p_b) & \text{if } D < 0 \end{cases}

DD' 的取值范围为 [1,1][-1, 1]D=1|D'| = 1 表示完全连锁不平衡,即四个单倍型中至少一个的频率为 0。

r2r^2 是 GWAS 中最常用的 LD 度量,定义为:

r2=D2pApapBpbr^2 = \frac{D^2}{p_A p_a p_B p_b}

r2r^2 的生物学直觉是:一个位点能在多大程度上预测另一个位点。具体来说,r2r^2 等于用一个位点的基因型预测另一个位点基因型时的决定系数。

r2r^2DD' 的关键区别:

  • D=1D' = 1 只要求四个单倍型中有一个缺失,但剩余三个单倍型的频率可能非常不均匀。
  • r2=1r^2 = 1 要求两个位点提供完全相同的信息量(即一个位点的基因型可以无误差地推断另一个位点)。

因此,r2r^2 对 GWAS 更有意义:它直接告诉我们在 GWAS 中检测到一个位点时,另一个位点是否也被”标记”了。

LD 最根本的来源是物理距离。两个位点在染色体上离得越近,减数分裂时它们之间发生重组(交叉互换)的概率就越低,关联就维持得越久。

数学描述:假设重组率为 cc(即每代两个位点间发生重组的概率),则 tt 代后 LD 系数的衰减为:

Dt=D0(1c)tD_t = D_0 (1-c)^t

cc 很小时(位点距离很近),(1c)tect(1-c)^t \approx e^{-ct},LD 以指数速度衰减。衰减到初始值一半所需的世代数为:

t1/2=ln2ct_{1/2} = \frac{\ln 2}{c}

对于人类基因组,c1 cM/Mbc \approx 1 \text{ cM/Mb}(即每兆碱基约 1% 的重组率),因此 1 cM 距离的两个位点需要约 69 代(约 2000 年)才能使 LD 衰减一半。

  • 瓶颈效应(Bottleneck):群体急剧缩小时,随机漂变会显著增强 LD。这是因为少数幸存者的染色体组合被”冻结”,随后的群体扩张使得这些组合在群体中广泛传播。
  • 选择(Selection):正向选择会带着周围的碱基一起扩散(选择性清扫 Selective Sweep),形成长距离的高 LD 区域。
  • 奠基者效应(Founder Effect):当一个小群体从大群体中分离出来建立新群体时,奠基者携带的有限单倍型多样性会产生广泛的 LD。

即使没有物理连锁,遗传漂变也可以在有限群体中产生 LD。设群体有效大小为 NeN_e,则在漂变作用下,DD 的期望衰减速度为:

E[Dt]=D0(112Ne)tE[D_t] = D_0 \left(1 - \frac{1}{2N_e}\right)^t

在小型群体中(NeN_e 较小),漂变产生的 LD 可以维持很长时间。这也是为什么人类群体中,非洲人群(历史上有效群体大小较大)的 LD 衰减比非非洲人群更快。

人类基因组的 LD 并不是平滑衰减的,而是呈现出块状结构

  • Block 内部:高度关联,仅需少数几个 Tag SNPs 就能代表整个区域的变异信息。
  • Block 边界:通常是重组热点(Recombination Hotspots),关联在这里迅速瓦解。

算法应用:这一特性是全基因组分型芯片设计的基础——我们不需要测序 30 亿碱基,只需检测几十万个精心挑选的 Tag SNPs 即可覆盖全基因组。

不同人群的 LD 衰减速度差异显著:

  • 非洲人群:LD 衰减最快,在约 10—20 kb 的距离上 r2r^2 即降至 0.2 以下。
  • 欧洲人群:LD 衰减中等,在约 20—50 kb 的距离上 r2r^2 降至 0.2。
  • 东亚人群:LD 衰减最慢,高 LD 可延伸至 50—100 kb 以上。

这种差异的历史原因是”走出非洲”事件:非洲以外的人群经历了奠基者效应和瓶颈效应,导致遗传多样性降低、LD 增强。

Tag SNP 的选择目标是找到最少数量的 SNP,使得每个非 Tag SNP 至少与一个 Tag SNP 的 r2r^2 超过阈值(如 r2>0.8r^2 > 0.8)。

这一问题是经典的集合覆盖问题(Set Cover Problem),已被证明是 NP-hard 的。实际工具采用贪心算法近似求解:

  1. 计算所有 SNP 对之间的 r2r^2 矩阵。
  2. 选择能”覆盖”最多未覆盖 SNP 的位点作为 Tag SNP。
  3. 重复直到所有 SNP 被覆盖。

基因型插补是利用 LD 关系,从已检测的位点推断未检测位点的基因型。

如果位点 A(已检测)与位点 B(未检测)之间的 r2=0.9r^2 = 0.9,那么我们可以用位点 A 的基因型来高置信度地推断位点 B 的基因型。

插补的统计模型通常基于隐马尔可夫模型(HMM),将参考面板(Reference Panel) 中的单倍型作为隐状态,被插补个体的基因型作为观测值。代表性工具包括 IMPUTE2、Minimac 和 Beagle。

插补的质量高度依赖于参考面板的规模和祖源匹配程度:

HapMap
早期的参考面板,包含约 1,400 个个体。规模较小,已被后续面板取代。
1000 Genomes Project
包含 2,504 个个体,覆盖 26 个种群。是目前最广泛使用的参考面板之一。
HRC (Haplotype Reference Consortium)
包含超过 32,000 个个体(主要为欧洲祖源)。由于样本量大,对欧洲祖源的插补质量远优于 1000 Genomes。
TOPMed
基于全基因组测序数据,包含超过 97,000 个个体。覆盖更多样化的祖源群体,是目前规模最大的参考面板之一。

GWAS 发现的显著信号通常覆盖一个 LD 区块中的多个 SNP。精细定位的目标是从这些连锁的 SNP 中识别真正的致病变异。

给定一个 GWAS 信号区域中的 mm 个 SNP,精细定位方法(如 CAVIAR、FINEMAP、SuSiE)计算每个 SNP 是因果变异的后验概率。核心思想是:

  • 利用 LD 矩阵对关联信号进行分解。
  • 在统计模型中允许存在多个独立的因果变异。
  • 通过后验包含概率(Posterior Inclusion Probability, PIP) 对 SNP 进行排序。

将 LD 信息与功能基因组注释相结合:

  • 位于启动子、增强子或转录因子结合位点中的 SNP 更可能是因果变异。
  • eQTL 数据可以帮助识别哪些 GWAS SNP 影响基因表达。
  • 保守性评分(如 PhyloP、GERP)可以评估变异的进化约束程度。

两个 SNP 位点,观察到以下单倍型频率:

单倍型频率
ABAB0.40
AbAb0.10
aBaB0.20
abab0.30

计算 DDDD'r2r^2

步骤 1:计算等位基因频率。

pA=0.40+0.10=0.50,pa=0.50p_A = 0.40 + 0.10 = 0.50, \quad p_a = 0.50 pB=0.40+0.20=0.60,pb=0.40p_B = 0.40 + 0.20 = 0.60, \quad p_b = 0.40

步骤 2:计算 DD

D=PABpApB=0.400.50×0.60=0.400.30=0.10D = P_{AB} - p_A p_B = 0.40 - 0.50 \times 0.60 = 0.40 - 0.30 = 0.10

步骤 3:计算 DD'

D>0D > 0,因此:

Dmax=min(pApb,papB)=min(0.50×0.40,0.50×0.60)=min(0.20,0.30)=0.20D_{\max} = \min(p_A p_b, p_a p_B) = \min(0.50 \times 0.40, 0.50 \times 0.60) = \min(0.20, 0.30) = 0.20

D=0.100.20=0.50D' = \frac{0.10}{0.20} = 0.50

步骤 4:计算 r2r^2

r2=D2pApapBpb=0.1020.50×0.50×0.60×0.40=0.010.060.167r^2 = \frac{D^2}{p_A p_a p_B p_b} = \frac{0.10^2}{0.50 \times 0.50 \times 0.60 \times 0.40} = \frac{0.01}{0.06} \approx 0.167

解读D=0.50D' = 0.50 表明两个位点之间存在中等程度的 LD(但不是完全连锁)。r2=0.167r^2 = 0.167 表明用一个位点预测另一个位点时,只能解释约 16.7% 的变异——这两个位点不能互相替代用于 GWAS 标记。