【基于R语言群体遗传学】-11-二倍体适应性进化

之前的博客，我们讨论了正选择在单倍体中的情况，可以先学习之前的博客：

群体遗传学_tRNA做科研的博客-CSDN博客

分子适应的经典观点是正选择。对于一个新的突变在进化意义上起作用，它必须在罕见时通过选择增加。之前例子中有选择系数s是负数，所以该等位基因有相对适应性损失，并且通常很快从种群中移除。这种类型的选择被称为负选择，虽然它乍看起来可能微不足道（迅速丢失的罕见突变对种群的一般或长期影响很小），但它是选择的一个非常重要的类别。

现在，让我们解决二倍体选择的问题。我们之前所有的选择示例都是基因选择：选择仅作用于单个等位基因。但是，由于二倍体有两个等位基因，它可以产生三种可能的基因型（AA、Aa、aa），每种基因型根据它们自己的选择系数可能有适应性优势或劣势：

sAA、sAa和saa

我们使用w = 1 + s来表示这三种基因型的相对适应性：

wAA、wAa和waa

因为我们处理的是基因型，我们可以通过将预期的基因型频率（根据哈代-温伯格方程）乘以它们的适应性来预测下一代中仅A等位基因的频率（p(t+1)）。我们必须记住，杂合子只对我们感兴趣的等位基因贡献一半的频率，因此我们需要将预期的杂合子频率2p(1−p)除以2以得到p(1−p)。

现在我们可以引入优势的复杂性。我们在之前讨论血型时简要提到了这个概念；具体来说，当A和B等位基因与i等位基因在杂合子中配对时，它们仍然分别赋予A型和B型血。所以我们可以说A和B等位基因对i等位基因是显性的，而i等位基因对A和B等位基因都是隐性的。就像表型一样，我们可以认为适应性效应是彼此显性和隐性的。

假设A等位基因有一个10%的适应性优势。我们的纯合子适应性值相当简单：

wAA = 1.1和waa = 1

如果A等位基因的适应性优势是显性的，我们的杂合子适应性将与AA纯合子相同：

wAa = wAA = 1.1

如果是隐性的，则与aa纯合子相同：

wAa = waa = 1

然而，显性/隐性关系并不一定如此明确。你可能会观察到一种称为不完全显性的现象，其中杂合子的适应性介于两个纯合子之间。在不完全显性的情况下，杂合子适应性是两个纯合子基因型适应性的确切平均值，我们将使用半显性（semidominance）这个词。这些显性程度预测了不同的进化轨迹。我们将使用相同的代码为所有三种适应性场景（隐性、显性和半显性）生成等位基因频率预测：

# 设置初始等位基因频率为 0.05
init_p <- 0.05

# 定义模拟的代数为 400 代
gen <- 400

# 定义三种不同的适应度场景：
# 隐性基因型的适应度值，纯合隐性基因型适应度为 1.1，杂合子适应度为 1
rec <- c(1.1, 1) 
# 显性基因型的适应度值，纯合和杂合适应度均为 1.1
dom <- c(1.1, 1.1) 
# 半显性基因型的适应度值，纯合子适应度为 1.1，杂合子适应度为 1.05
sem <- c(1.1, 1.05) 

# 创建一个矩阵，包含三个相同的初始等位基因频率
p <- matrix(c(init_p, init_p, init_p))

# 将三种适应度场景组合成一个列表
w <- list(rec, dom, sem)

# 定义一个函数 FitFreq，用于计算下一代的等位基因频率
FitFreq <- function(X, p){
  # 计算总适应度
  w_total <- X[1]*p^2 + X[2]*2*p*(1-p) + (1-p)^2
  # 根据适应度计算下一代的等位基因频率
  p_t <- (X[1]*p^2 + X[2]*p*(1-p)) / w_total
  # 返回计算出的下一代等位基因频率
  return(p_t)
}

# 使用隐性基因型的适应度值和初始等位基因频率调用 FitFreq 函数
# 以计算在这些条件下下一代的等位基因频率
FitFreq(w[[1]], p[[1]])

我们得到：

我们对三种情况进行建模模拟：

# 循环从第1代到第399代（总共400代，减去初始的1代）
for(i in 1:(gen-1)){
  # 使用cbind函数将新计算的等位基因频率添加到矩阵p中
  # lapply函数遍历适应度场景列表w，并对每个场景应用FitFreq函数
  # seq_along(w)生成与w长度相同的序列，作为lapply的索引
  # function(j, y, n) {FitFreq(y[[j]], n[[j]])}是一个匿名函数，
  # 它接受索引j，并使用w和p的最新列来计算新的等位基因频率
  p <- cbind(p, lapply(seq_along(w),
                       function(j, y, n) {FitFreq(y[[j]], n[[j]])},
                       y=w, n=p[,ncol(p)]))
}

# 绘制图形，设置x轴标签为"Generations"，y轴标签为"Allele frequency"
# x轴范围设置为1到400代，y轴范围设置为0到1
plot(x=NULL, xlab="Generations", ylab="Allele frequency",
     xlim=c(1,gen), ylim=c(0, 1))

# 定义颜色和线型
colors <- c("orange", "black", "cyan")
line <- c(1,2,4)

# 循环遍历矩阵p的每一行（每种适应度场景）
for(i in 1:nrow(p)){
  # 绘制线条，表示每一代中该适应度场景下的等位基因频率变化
  # lwd设置线宽，lty设置线型，col设置颜色
  lines(1:gen, p[i,], lwd=2, lty=line[i], col=colors[i])
}

# 添加图例，说明每条线代表的适应度场景
legend("bottomleft",
        legend=c("Recessive","Dominant","Semi-dominant"),
        inset=c(0,1), xpd=TRUE, bty="n",
        col=colors, lty = line, lwd=2)

我们可以看到，当基因频率罕见时，显性适应度效应在频率上的初始上升速度更快，但在高频率下固定的速度较慢；隐性适应度效应则相反，在很长一段时间内频率几乎没有增加，然后迅速上升到固定。半显性适应度效应似乎是在整体固定时间方面最有效的二倍体配置。然而，直接作用于单倍体基因型的基因选择仍然是最有效的选择形式。

半显性Logit预测

如果适应度效应是半显性的，那么通过logit变换（即取几率的自然对数）可以线性化预测的等位基因频率变化。在广泛的频率范围内，其他类型的显性在通过logit变换后几乎是线性的。

这使我们能够使用线性回归从等位基因频率数据中估计选择系数。让我们看看Fisher和Ford（1947）的真实时间等位基因频率测量结果，这些结果显示了猩红虎蛾（Callimorpha dominula）八代中等位基因频率的变化。

所讨论的等位基因（在89-90%的频率下相当常见）影响蛾的翅膀色素沉着，并可能由于捕食（花纹可能会暴露你）或配偶偏好和繁殖成功（某些颜色图案比其他颜色更具吸引力）或你能想到的任何其他选择故事而影响生存。我们可以看到常见的纯合子的翅膀是红色和黑色，带有黄色斑点，而较罕见的纯合子的斑点要少得多，翅膀通常更暗。此外，我们实际上可以看到，我们所测量的等位基因是不完全显性的，因此杂合子与两种纯合子都有区别，具有中间表型。

# 加载 popgenr 包
library(popgenr)

# 载入蛾类种群数据
data(moth)

# 从数据中提取年份和等位基因频率
Years <- moth$Years
Freq <- moth$p

# 计算等位基因频率的对数几率（logit 转换）
logit <- log(Freq/(1-Freq))

# 绘制年份与对数几率的关系图
plot(Years, logit)

# 进行线性回归分析，拟合年份与对数几率之间的关系
(linear <- lm(logit ~ Years))

# 在图中添加线性回归线
abline(linear)

很明显等位基因的频率随着时间的推移而增加。回归最佳拟合线的斜率实际上被假定为大约等于携带一个等位基因副本的杂合子和携带两个副本的纯合子之间的适应度差异。在我们的输出中，我们看到斜率约为0.12；这为我们提供了杂合子的选择系数的估计值（s+− ≈ 0.12），如果我们假设半显性，那么我们就可以估计我们剩余的基因型适应度为w++ = 1.24，w+− = 1.12，和w−− = 1，其中“+”表示我们在数据中追踪的等位基因，“-”是导致较暗表型的较罕见等位基因。但这种方法本身并不能区分自然选择和遗传漂变。我们假设选择是我们等位基因频率变化的主要驱动力，但完全有可能漂变也在起作用。此外，这些测量并不是相互独立的——明年的等位基因频率取决于前一年的频率

一个非常灵活的测试数据集框架是在模型下生成模拟数据，在这个例子中是一个只有遗传漂变而没有选择的零模型，并观察观测结果在模拟结果的分布中的位置。让我们首先保存来自线性回归的平方相关系数，并将其作为描述我们的等位基因频率随时间变化的关键描述:

(obs_r <- summary(linear)$r.squared )
[1] 0.6551432

我们可以在漂变模型下模拟数千个数据集，并计算模拟数据集超过这个值的频率。我们甚至可以允许我们的种群大小和选择系数都发生变化，以估计s+−的可能值范围；修改之前的遗传漂变代码，我们可以看到，如果实际上根本没有选择，只有漂变在起作用，我们有多频繁地看到相同或更大的r方值。种群大小在这里将发挥重要作用；我们已经看到，随着我们的种群规模变小，漂变的影响变得更加明显，并且很容易淹没选择的影响。因此，如果我们用小种群规模运行这个模拟，我们将放大漂变的影响。为了防止这种情况，我们可以说我们正在研究的种群至少有1000个个体。

# 设置初始等位基因频率
initp <- 0.908

# 定义模拟的代数
gen <- 8

# 创建一个向量，表示年份
Years <- c(1:gen)

# 设置模拟重复的次数
reps <- 10000

# 初始化计数器
count <- 0

# 设置种群大小
N <- 1000

# 开始模拟
for(i in 1:reps){
  # 重置等位基因频率
  p <- initp
  # 进行除第一代外的每一代的模拟
  for(j in 1:(gen-1)){
    # 使用二项分布模拟等位基因在下一代中的传递
    a <- rbinom(1, 2*N, p[j])
    # 更新等位基因频率
    p <- c(p, a/(2*N))
  }
  # 对等位基因频率进行对数几率转换
  logit <- log(p/(1-p))
  # 拟合线性回归模型
  reg <- lm(logit ~ Years)
  # 获取模拟得到的 r^2 值
  sim_r <- summary(reg)$r.squared
  # 如果模拟得到的 r^2 值大于等于观测到的 r^2 值（obs_r），则计数器加一
  if(obs_r <= sim_r){
    count = count + 1
  }
}

# 计算满足条件的模拟次数占总模拟次数的比例
count / reps

从这个结果中，我们可以看到，由于漂变看到观测数据的概率大约是33%。这个比例足够高，以至于我们应该合理地对关于这个等位基因的任何选择故事持怀疑态度。在我们继续讨论其他形式的选择之前，关于这个数据集的最后一点是，早期采集的猩红虎蛾样本中，常见等位基因的频率甚至更高。所以，在某个时间点，人们实际上认为选择是在增加“-”等位基因的频率，而不是像我们在1939年到1946年看到的减少它。这应该强调在试图推断选择时应谨慎行事，因为除了用随机漂变解释变化外，选择压力（以及种群大小）也可能随时间波动。

下一篇文章我们将考虑一种有意思的现象-杂合优势。