连锁不平衡的计算以及LDSC分析多基因遗传

连锁不平衡（linkage disequilibrium,LD）是指在某一个群体中，不同座位上两个基因同时遗传的频率明显高于预期的随机频率现象，连锁不平衡的程度通常用 r2 来衡量。

D是LD的基本单位，度量观察到的单倍型频率与平衡状态下期望频率的计算方法如下：

D=P(AB)-P(A)*P(B)

P(AB)表示实际观察到的AB频率，P(A)*P(B)表示AB频率的期望值。（如果发生连锁不平衡，实际观测到的AB频率肯定不等于AB频率的期望值）

如果D值显著偏离0，则说明存在LD。因为D的取值强烈地依赖于人为制定的等位基因频率，所以它不利于LD程度的比较。标准化的不平衡系数D'能够避免这种对

等位基因频率的依赖。D'的计算方法如下：

D'=D/Dmax

当D<0, Dmax=min{P(A)P(B),P(a)P(b)}；

当D>0, Dmax=min{P(A)P(b),P(a)P(B)}；

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

当D‘=0，表示连锁完全平衡，随机组合；

除了D值之外，还有一个衡量连锁不平衡程度的标准，就是`r2`值，计算公式如下

r2=D*D/(P(A)P(a)P(B)P(b))

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

当r2=0，表示连锁完全平衡，随机组合

下面通过一个示例，看下实际分析中，P(A)和P(B)如何计算。

在一个群体中，观测到的单倍型频率分布如下

根据单倍型的频率分布，可以计算出如下的等位基因频率

单倍型A1B1的D值计算如下

对于不同单倍型，其D值汇总如下

上面的例子，可以帮助我们理解D值计算的过程。在实际分析中，我们通常会拿到样本的基因分型文件，通过这个文件我们可以非常容易的计算出allel的频率，但是对于单倍型的频率是不能直接计算得到的，都是借助算法估算出单倍型的概率，然后进行计算。对于连锁不平衡的计算，有很多的软件可以用，最常用的有plink和haploview，当然也有很多的R包可以进行计算。

理论上来说任意两个基因之间都可能存在连锁不平衡，但是实际操作中，认为只有一定区间范围内的基因会存在连锁不平衡，距离大于区间的基因，两者出现连锁不平衡的概率非常小，所以就不去计算。对于连锁不平衡的结果，通常采用heatmap热图的形式进行展示，haploview给出的LD heatmap 示例如下

颜色从白色到红色，代表连锁程度从低到高，方框中的数值为r2，为了美观，这里将r2 乘以了100。除了heatmap 之外，还会有下面这种散点图

横坐标为基因之间的距离，纵坐标为衡量连锁不平衡的R2值。从图中可以看出连锁不平衡的规律，在一定距离内存在连锁不平衡程度较高，大于一定距离后，出现连锁不平衡的概率就大大降低了。这就是为什么在实际操作中只计算一定范围内的连锁不平衡的原因。

LDSC（连锁不平衡回归分析）方法校正GWAS分析中的多基因遗传效应

GWAS研究中，混杂因素偏倚（隐晦的亲缘关系和种群分层）和多基因效应是造成检验统计量膨胀的主要原因。

通过引入LDSC（ Linkage Disequilibrium Score Regression）检测统计值与LD（r2）之间的回归关系，可以量化每个变量的贡献，回归分析的截距值也可以用来区分统计量膨胀的原因并对进行校正。

引入如下模型，变异 j 的预期的卡方检验统计量表示如下：

N：样本大小，M：SNP数量，h2/M：每个SNP的平均遗传率（遗传率的计算方法下一个笔记讲），a：混杂偏倚的效应值

由此，回归方程中的截距值-1就是混杂偏倚效应对卡方统计量膨胀贡献均值的估计值。

利用LDSC工具绘图，若得到左图（截距值近似于1），表示不存在混杂效应偏倚；若得到右图，表示可能存在混杂效应偏倚。

参考文献：https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4495769/

这篇文献还没看完全看懂，等我来补QAQ