[连锁不平衡粗俗的说就是：这几个基因耍流氓，喜欢抱团遗传，不再随机。而连锁不平衡衰减是指在基因组上，随着物理距离的增大，两个连锁的的等位基因的连锁程度不断减小。]

LD衰减 图，在重测序类的文章中会经常出现群体遗传、GWAS等的文章里面 。

要理解LD衰减图，我们就必须先理解连锁不平衡（Linkage disequilibrium，LD）的概念。 连锁不平衡是由两个名词构成，连锁+不平衡。 前者，很容易让我们产生概念混淆；后者，让这个概念变得愈加晦涩。因此从一个类似的概念入手，大家可能更容易理解LD的概念，那就是基因的共表达。

基因的共表达，通常指的是两个基因的表达量呈现相关性。 比较常见的例子就是：转录组因子和靶基因间的关系。因为转录因子对它的靶基因有正调控作用，所以转录因子的表达量提高会导致靶基因的表达量也上调，两者往往存在正相关关系。这个正相关关系，可以使用相关系数r^2来度量，这个数值在-1~1之间。总而言之，相关性可以理解为两个元素共同变化，步调一致。

类似的，连锁不平衡（LD）就是度量两个分子标记的基因型变化是否步调一致，存在 相关性 的指标。如果两个SNP标记位置相邻，那么在群体中也会呈现基因型步调一致的情况。比如有两个基因座，分别对应A/a和B/b两种等位基因。如果两个基因座是相关的，我们将会看到某些基因型往往共同遗传，即某些单倍型的频率会高于期望值。

例如在下图2中，在群体中（A，a，B，b）各个基因型的频率已知的情况下，各种单倍型的期望频率（AB、Ab、aB、ab）都是可以计算出来。例如，AB的频率=（A的频率）X（B的频率）。但我们实际统计群体中各个单倍型的频率的时候，会观察到某些单倍型的频率会大于期望值，例如下图中的单倍型AB的理论频率是0.12，但观察到的实际频率是0.29。那么说明，基因型A更倾向于基因型B共同遗传。

这一般往往是由于在祖先的基因组中，A和B就是位于同一条染色体上，在传代过程中，这种共同遗传的关系被保留了下来。位点间的这种相关性，在杂交家系中一般被称为连锁（孟德尔老师豌豆实验中的发现），在自然群体中则一般被称为连锁不平衡。所以连锁不平衡中的“不平衡”，我认为可以理解为单倍型的频率分布偏离期望值，偏离了平衡。

这种不同基因座间的相关性，用一个数值来衡量就是D值（图2中有计算公式）。类似相关系数是标准化后的协方差，LD系数（r 2）则是标准化后的D值（图2中有计算公式），这个数值在0~1波动。r 2=0就是两个位点完全不相关，群体中单倍型分布是随机的（观测值=期望值）。r^2=1就是两个位点完全相关，某些基因型（A）只与特定的基因型（B）共同出现。

一般而言，两个位点在基因组上离得越近，相关性就越强，LD系数就越大。反之，LD系数越小。也就是说，随着位点间的距离不断增加，LD系数通常情况下会慢慢下降。这个规律，通常就会使用LD衰减图来呈现。

图形的解读

LD衰减图就是利用曲线图来呈现基因组上分子标记间的平均LD系数随着标记间距离增加而降低的过程。大概的计算原理就是先统计基因组上两两标记间的LD系数大小，再按照标记间的距离对LD系数进行分类，最终可以计算出一定距离的分子标记间的平均LD系数大小。如图3是黄瓜重测序文章中统计各个亚群体的LD衰减速度的图形。横坐标是物理距离（kb），纵坐标是LD系数（r^2）。

从图中我们可以看出，西双版纳这个亚群体（紫色线）在基因组上50kb距离的平均LD系数大小约为0.4，但到了100kb的距离，对应的平均LD系数大小则降低到了不到0.3。而且，我们从图中也可以观察到LD系数的衰减速度在不同的亚群体快慢不同，衰减速度是 india >East Asian&Eurasian >Xishuanbanna。那说明india群体的LD衰减距离最小，可能是india这个群体遗传多样性最高导致。这句话该如何理解呢？

实际上，LD衰减的速度在不同物种间或同物种的不同亚群体间，往往差异非常巨大。所以，通常会使用1个标准——“LD衰减距离”来描述LD衰减速度的快慢。

LD衰减距离通常指的是：当平均LD系数衰减到一定大小的时候，对应的物理距离。

“一定大小”是这个定义的关键点，但没有特别统一的标准，在不同文章中标准不同。常见的标准包括：a）LD系数降低到最大值的一半；b）LD系数降低到0.5以下；c）LD系数降低到0.1以下；d）LD系数降低到基线水平（但注意，不同材料的基线值是不同的。比如图3黄瓜群体的基线大概是0.1）。

所以，下次你在文章中看到“LDdecay distance is XXkb”的时候，别忘了看看作者使用的标准是什么。

如图3所示， LD系数衰退速度会受到不同因素的影响而有所不同。常见的因素包括：

1）物种类型LD存在的本质是两个位点的连锁遗传导致的相关性。 但这种相关性理论上会随着世代的增加、重组次数的增加而不断下降。所以，那些繁殖力强、时代间隔短的物种（例如，昆虫），其LD衰减的速度是非常快的。例如在家蚕和野蚕群体中，LD系数下降到最大值的1/2仅仅需要46bp和7bp的距离[3]。

2）群体类型相同物种的不同群体，由于其遗传背景不同，LD衰减速度也存在很大的差异 。驯化选择，会导致群体遗传多样性下降，位点间的相关性（连锁程度）加强。所以，通常驯化程度越高，选择强度越大的群体，LD衰减速度是最慢的。例如，栽培稻比野生稻通常更大的LD衰减距离。类似的，自然选择、遗传漂变导致的群体遗传多样性下降，也会减慢LD衰减的速度。

3）在染色体的位置染色体不同区域的LD衰减距离而是不同的。 通常着丝粒区更难重组，所以LD衰减更慢。而基因组上那些受选择的区域相比普通的区域，LD衰减速度也是更慢的[3]。

LD衰减速度，在群体遗传分析中本身是对群体特性的评估，与群体类型的特性（自然群体还是驯化群体，选择强度大小）是相关的。但在其他研究中还有更多的应用价值。

基于分子标记（例如，SNP芯片，GBS测序）的GWAS分析，其实并没有检测到功能突变，本质就是利用标记和功能突变的相关性（LD关系），来检测与性状相关的功能突变的位置。一般而言，LD系数大于0.8就是强相关。如果LD系数小于0.1，则可以认为没有相关性。如果LD衰减到0.1这么大的区间内都没有标记覆盖的话，即使这个区间有一个效应很强的功能突变，也是检测不到关联信号的。所以，通常可以通过比较LD衰减（到0.1）距离和标记间的平均距离，来判断标记是否对全基因组有足够的覆盖度。

而如果GWAS检测到显著关联的区间后，则可以通过进一步绘制局部的LD单体型块图，来进一步判断显著相关的SNP和目标基因间是否存在强LD关系。这个图形我们下一篇文章会介绍。

再提一个应用的例子。在之前的文章中我们提到过，在进行STRUCTURE分析的时候理论上必须输入不相关的位点。那么，就可以通过预估LD衰减到0.1的距离，来判断标记间的距离必须大于多少才能保证标记间不具相关性（LD<0.1）。

3.绘制方法

LD衰减图的绘制，实际上有两个步骤：

1）计算marker间两两的LD系数大小

这个可以使用haploview软件完成。计算的时候，只要设定一个关键的参数：区间大小。例如设定为5Mb，那么软件就会计算基因组上所有距离<5Mb的两两位点间的LD系数。实际上这个参数设定更大也没有意义，一般情况下位点间的相关性不会延伸到大于5Mb这么远的距离。

2）绘图

将LD系数按照对应的两个marker间的距离进行分类，例如：距离按照区间大小0 5k,5k 10k,10k~15k…..分别分类。如果重测序的数据，SNP标记密度较大，这个分类区间可以设置小一些；如果是简化基因组数据，SNP标记较为稀疏，则分类区间可以适当加大。然后计算每种距离分类的LD系数的均值。最后在利用均值绘制曲线图就ok了。这一步的绘图，使用excel或R语言都可以轻松完成。

连锁不平衡 （linkage disequilibrium, LD）分析是群体遗传学研究中常见的分析内容，也是关联分析的基础，在很多的GWAS文章中都会出现LD衰减图及单倍型block图，接下来一起连锁不平衡（linkage disequilibrium, LD）初探。

图1 水稻自然群体连锁不平衡衰减与桃果糖含量位点相关单倍型区块

1、 LD的概念

当位于某一座位的特定等位基因与另一座位的某一等位基因同时出现的概率大于群体中因随机分布的两个等位基因同时出现的概率时，称这两个座位处于连锁不平衡状态。

2、 LD的计算方法与度量指标

2.1 D值的计算

LD的基本单位是D值，度量观察到的单倍型频率与平衡状态下期望频率的偏差。D值根据单倍型频率必≥0，计算取值范围为[-0.25,0.25]。

D=Pr（A，B）-Pr（A）×Pr（B）

=PAB-PAPB

=PAB-（1-Pa）（1-Pb）

=PAB-（1-Pa-Pb+PaPb）

=PAB-（PA-Pb+PaPb）

=PAB-[PAB+PAb-（PAb+Pab）+PaPb）]

=PAB-（PAB-Pab+PaPb）

=Pab-PaPb

2.2 标准化指标：D’和r2

由于D值严格依赖于等位基因频率（allele frequency），故不适合应用于表述实际的LD强度，最常用度量LD的是D’和r2，两者都基于D。D’反应群体的重组历史，适用于研究群体连锁不平衡程度，r2反应等位基因相关程度，适用于GWAS。LD衰减作图中通常采用r2来表示群体的LD水平；Haplotype Block中通常采用D'来定义Block；迁移、突变、选择、有限的群体大小以及其他引起等位基因频率改变的因素，这些都会引起LD变化，下面尝试计算一下LD吧~

图2 D’和r2的计算

小Tips ：

当D'=0，r2=0时，处于完全连锁平衡状态

当D'=1，r2=1时，处于完全连锁不平衡状态。

其中，从0—1之间的度量越高，LD越高，如果两个位点连锁，连锁程度也越高。

3、 LD衰减

LD衰减指位点间由连锁不平衡到连锁平衡的演变过程；LD衰减的速度在不同物种间或同物种的不同亚群间，差异非常大。所以，通常会使用“LD衰减距离”来描述LD衰减速度的快慢，不同文章中“LD衰减距离”标准不同，常见的标准有：LD系数降低到最大值的一半、LD系数降低到0.5以下、LD系数降低到不同物种的基线水平等，我们在阅读文献时有必要留意文章使用的标准~

LD衰减距离在群体遗传学中的应用也非常广泛，一方面可以判断GWAS所需标记量，决定GWAS的检测效力以及精度；另外也可以辅助分析进化与选择，在同一个连锁群上，LD衰减慢说明该群体受到选择，一般来说，野生群体比驯化改良群体LD衰减快，异花授粉植物比自花授粉植物LD衰减快。

图3 LD下降到最大值一半对应的物理距离

4、 LD分析的软件与画图命令

目前比较常用的计算 LD 的软件为 PLINK ，但不支持直接读取 VCF 格式的文件，使用 PLINK 计算 LD 之前需要先将 VCF 格式的文件转换为 PED 格式或 bed + bim + fam 的格式。这样的格式转换会造成额外的存储负担。而另一款软件PopLDdecay 一个主要的优点在于可以读取 VCF 格式的文件，直接生成 LD 统计数据并画出 LD 的衰减图。以软件PopLDdecay计算绘制LD图的命令如下：

#PopLDdecay -InVCF Final_snps.vcf -OutStat Out.LDdecay -SubPop pop.list

全基因组关联分析（Genome-Wide Association Study， GWAS）是以连锁不平衡（LD）为基础，利用全基因组范围内群体中高密度的分子标记，鉴定与复杂性状表型变异相关联的分子标记，进而挖掘与表型相关基因的方法。

D' 的计算方法如下：

当 D <0, D max = max { - P A P B , -(1-P A ) (1- P B ）} ；

当 D >0, D max = min { P A (1- P B ）, (1- P A ) P B } ；

当 D ' =1，表示连锁完全不平衡，没有重组；

当 D ' =0，表示连锁完全平衡，随机组合，独立遗传。

但 D ' 也有它的局限性，比如当单倍型为两种或三种时，| D ' |一定等于1，但是当| D ' |<1时， D ' 的值究竟表示多大程度的连锁不平衡，是很难做出准确判断的。另外 D ' 严格依赖于样品的大小，如果样本偏少时，SNP数量比较少，这样算出来的 D ' 就会偏大，尤其是某个位点其中一个等位基因频率很低时，因此较高 D ' 背后，实际上可能是连锁不平衡程度很低的两个位点。因此引进 r 2 来表示LD， r 2 的计算方法如下：

当 r 2 =1，表示连锁完全不平衡，没有重组；

当 r 2 =0，表示连锁完全平衡，随机组合，独立遗传。

---