【基于R语言群体遗传学】-16-中性检验Tajima‘s D及连锁不平衡 linkage disequilibrium (LD)

Tajima's D Test

已经开发了几种中性检验，用于识别模型假设的潜在偏差。在这里，我们将说明一种有影响力的中性检验，即Tajima's D（Tajima 1989）。Tajima's D通过比较数据集中的两个𝜃 = 4N𝜇估计值来工作。我们已经推导出了𝜃is，等于平均成对杂合性（average pairwise heterozygosity），当我们讨论共祖时（也称为Tajima的估计器）。当考虑DNA序列集合中的等位基因或SNP总数以及将它们联合到共同祖先的共祖树内包含的预期世代数时，可以得出另一种推导。这被称为Watterson的𝜃估计器或𝜃W（Yong 2019）。正如我们在共祖中所展示的，n个谱系的预期共祖时间是：

一组n个初始谱系到一个单一祖先的所有共祖时间的总和是：

在每次共祖事件之间的步骤中，有i + 1个谱系可能发生突变。因此，当考虑一段时间内所有谱系可能产生的等位基因数量时，我们乘以i + 1：

我们将树上所有谱系的总时间乘以每代突变率𝜇，得到我们期望在n个DNA序列样本中的等位基因总数，而4N是一个常数，所以我们可以将其放在求和之外

其中S是序列集合中SNP的数量。这可以重新排列以从SNP的数量估计𝜃W=4N𝜇：

请注意，Watterson的𝜃估计器需要了解一组谱系的共祖，但这是在1975年发表的，当时还没有发表超过两个谱系的共祖（Kingman 1982）。Tajima's D是平均成对杂合性𝜃估计值与从样本中SNP数量估计的𝜃之间的差异，除以该差异的预期方差的平方根：

而：

这看起来有点乱，但在各种位置反复出现两个不同的n求和，只需要计算一次然后填入。这可以用以下代码计算（平均成对差异和S也可以从数据集中计算，但为了简洁起见，我们在这里省略了）。 

Tajima's D是用来评估一个种群中中性突变（即没有自然选择影响的突变）的假设是否成立。它通过比较两个不同的估计器来衡量种群的遗传多样性和种群规模的变化，我们通过R语言实现：
 

# Calculates Tajima’s D# 平均成对差异，用于估计theta_IS
theta_IS <- 2.8# 数据集中的SNP数量
S <- 16# 采样的等位基因拷贝数
n <- 20# 初始化求和变量
i1_sum <- 0.0# 循环计算i1的和，这是Watterson's theta的一部分
for(i in 2:n-1){i1_sum <- i1_sum + 1/i
}# 计算Watterson's theta，它是基于序列多态性的一个种群规模的估计器
theta_W <- S / i1_sum# 初始化第二个求和变量
i2_sum <- 0.0# 循环计算i2的和，用于后续计算
for(i in 2:n-1){i2_sum <- i2_sum + 1/i^2
}# 计算期望值e1，它是Tajima's D公式中的项
e1 <- ((n+1)/(3*(n-1)) - 1/i1_sum) / i1_sum# 计算期望值e2，它也是Tajima's D公式中的项
e2 <- (2*(n^2+n+3)/(9*n*(n-1)) - (n+2)/(n*i1_sum) + i2_sum/i1_sum^2) / (i1_sum^2 + i2_sum)# 计算Tajima's D值，它衡量的是theta_IS和theta_W之间的标准化差异
(D <- (theta_IS - theta_W) / sqrt(e1*S + e2*S*(S-1)))

前三个变量将根据您的数据进行调整。在这个例子中，返回的D = -1.409。大于或小于2的D被认为是显著的；然而，实际的p值是通过模拟确定的。D的正值表示中间频率等位基因过多，这可能是由于人口减少或平衡选择，因为这两种情况都会延长人口历史较老部分的共祖事件时间。在更大的种群中，有更多的祖先可供选择，共祖是一种罕见的事件，并会膨胀𝜃IS，因为相对于𝜃W，较老的谱系在后代中以更高的频率共享。这种负D表示稀有频率等位基因过多（共祖的最近尖端被放大；它们对S和𝜃W的贡献比对平均成对差异的贡献更多，因为它们很稀有），这表明人口扩张、选择性清除或对有害等位基因的低效净化选择。关键是要计算多个位点的D值，并寻找异常值以标记假定选择候选者。人口统计学效应，如人口规模的变化，应该影响基因组中的所有位点，而选择（通常被认为）在其影响上是位点特异性的。还有许多其他的中性检验，如HKA检验（Hudson等人，1987）、McDonald-Kreitman检验（McDonald和Kreitman，1991）、Fay和Wu的H（Fay和Wu，2000）以及dN/dS比率（Yang和Bielawski，2000）。其中许多也利用了物种之间发生的遗传变化，它们都有各自的优点和缺点。

linkage disequilibrium (LD)

群体遗传学的一个独特性质，在进化博弈论等类似领域中并未发现，即不同位点甚至不同染色体上的等位基因可以“链接”（尽管并不总是与经典遗传学中的重组图谱同义），并且比随机预期更频繁地一起遗传。 连锁不平衡（LD）的程度由𝒟量化，不要与Tajima的D混淆。考虑两个位点：一个具有A/a多态性，另一个具有B/b多态性。我们对于跨位点一起遗传的等位基因之间的关联感到好奇。使用概率的乘法规则，我们期望AB单倍型（pAB）的频率是两个等位基因频率pApB的乘积，如果它们是独立遗传的话。这两者之间的差异由𝒟量化，作为连锁不平衡的一种度量。

根据AB单倍型是过量还是不足，𝒟可以是正数或负数（或者如果你从ab单倍型任意计算𝒟，符号会改变）。𝒟也可以从所有单倍型频率计算得出。

为了说明，假设我们有一个包含两个SNP的单倍型频率的小型数据集。一个是A/G多态性，另一个具有C/T等位基因：

让我们关注A-C单倍型。A等位基因的频率是0.57，C等位基因的频率是0.65

遗传漂变、种群结构和强选择是推动𝒟偏离零的力量。在大种群中，预测𝒟随时间呈指数衰减轨迹返回到零，就像在小种群中遗传漂变下的杂合性一样：

其中r是感兴趣的位点对之间预期的重组分数。这可以用来估计单倍型的年龄。最后，即使对于在不同染色体上独立分配的位点，𝒟也需要时间衰减。哈代-温伯格基因型可以在一代中恢复，但过去的事件对LD有持续影响，这可以用来推断更远的过去的过程，如种群中不再存在的种群结构。 当从实际数据集计算𝒟时，双重杂合子是不明确的。假设我们有一个个体的C/T，A/G SNP集合。C等位基因与第二个位置的A还是G等位基因相关联？通常我们不知道。但是，不明确的单倍型频率为我们提供了关于解决双重杂合子可能方式的信息。如果C-G单倍型非常常见，而C-A单倍型很少见，那么这表明C/T，A/G个体可能具有C-G/T-A单倍型。使用这种方法计算𝒟太繁琐，无法手工完成。 幸运的是，这正是EM算法发挥作用的地方。Kalinowski和Hedrick（2001）使用大角羊（Ovis canadensis）数据集（Boyce等人，1997）来估计LD。这个物种很罕见，样本量很小，所以我们需要从可用的数据中获得尽可能多的信息。 以下R代码实现了Kalinowski和Hedrick（2001）给出的方程式。它从猜测相等的单倍型频率和𝒟 = 0开始。然后它更新这个猜测，并迅速达到最大似然解𝒟 ≈0.0779和na-B单倍型频率基本为零。

# 定义各个复合基因型的频率
AABB <- 2  # 两个位点都是纯合子AABB的个体数量
AaBB <- 0  # 一个位点是杂合子，另一个是纯合子AaBB的个体数量
aaBB <- 0  # 一个位点是纯合子aa，另一个是纯合子BB的个体数量
AABb <- 0  # 第一个位点是纯合子AA，第二个位点是杂合子Bb的个体数量
AaBb <- 1  # 两个位点都是杂合子AaBb的个体数量（双杂合子）
aaBb <- 0  # 第一个位点是纯合子aa，第二个位点是杂合子Bb的个体数量
AAbb <- 1  # 第一个位点是纯合子AA，第二个位点是纯合子bb的个体数量
Aabb <- 0  # 第一个位点是杂合子Aa，第二个位点是纯合子bb的个体数量
aabb <- 0  # 两个位点都是纯合子aabb的个体数量# 使用上述输入运行函数‘Dcalc’
Dcalc(AABB, AaBB, aaBB, AABb, AaBb, aaBb, AAbb, Aabb, aabb)

𝒟与用于衡量线性相关性的统计相关系数（皮尔逊）“r”有关。让我们使用ℛ表示相关系数，以避免将其与重组分数r混淆。𝒟2除以所有等位基因频率的乘积等于ℛ^2。

此外，奇怪的是，如果我们将ℛ2乘以采样的染色体总数（如果我们观察n个二倍体个体，通常为2n），那么我们就会得到一个具有一个自由度的𝜒2统计量：

然而，这并不令人惊讶。在如此小的样本量下，即使LD非常强，检测偏差的能力也非常有限。最后，我们希望指出，EM算法是一种“爬山”算法，它找到一个局部最大似然峰值。可能存在其他峰值，可以使用MCMC方法来处理这个问题，并更全面地探索复杂的似然表面。