【基于R语言群体遗传学】-8-代际及时间推移对于变异的影响

上一篇博客，我们学习了在非选择下，以二项分布模拟遗传漂变的过程：【基于R语言群体遗传学】-7-遗传变异（genetic variation）-CSDN博客

那么我们之前有在代际之间去模拟，那么我们就想知道，遗传变异随着时间的推移，会发生什么样的变化，我们今天进行探讨：

遗传漂变的基准实证研究之一是由早期遗传学对于果蝇的研究提供的。在五十年代，彼得·布里进行了一个令人惊讶的费力实验，让超过100个小型果蝇群体独立进化，这些果蝇具有突变眼色表型，持续了二十代。他以50%的频率开始计算每组16只果蝇，携带两个bw（“棕色”）等位基因，这使纯合子个体呈现出鲜艳的红橙色眼睛，而杂合子个体则可以通过明显较浅的橙色眼睛来区分。不携带bw75等位基因（bw/bw纯合子）的果蝇会表现出白色眼睛。为了保持种群大小不变，每一代都从上一代的后代中随机抽取一组新的16只果蝇，作为下一代的父母。通过观察眼睛颜色，布里统计了超过100个重复实验中bw75拷贝数随时间的变化，而且整个过程做了两次！

在我们查看布里的数据之前，让我们首先以一种新的方式可视化我们之前的模拟，以便我们可以更容易地进行比较，然后再使用我们模拟的数据和布里的数据可视化进行比较：

# 设置群体大小
N <- 16

# 可能的等位基因数量范围（从0到2倍群体大小）
possible <- 0:(2*N)

# 初始化一个矩阵x，用于存储概率值，初始时所有概率为0
x <- matrix(c(rep(0,2*N+1)), ncol=2*N+1, byrow=T)

# 在矩阵x的中心位置（即N+1的位置）设置概率为1，代表起始状态
x[N+1] <- 1

# 初始化一个空向量P，用于存储二项分布的概率
P <- NULL 

# 计算每个可能的等位基因数量的二项分布概率，并存储在P中
for(i in possible){
  P <- c(P,dbinom(possible, size=2*N, prob=i/(2*N)))
}

# 将P转换为矩阵Q，每行对应一个可能的等位基因数量的所有概率
Q <- matrix(P, ncol=2*N+1, byrow=T)

# 计算初始状态下的概率分布R
R <- x %*% Q 

# 将初始概率分布赋值给Prob
Prob <- R

# 初始化代数计数器g
g <- rep(1,ncol(R))

# 循环直到达到第19代
while(g[1]<19){
  # 更新概率分布R
  R <- R %*% Q
  # 将新的概率分布添加到Prob矩阵中
  Prob <- rbind(Prob,R)
  # 代数计数器加1
  g <- g+1
}  

# 使用persp函数绘制概率分布的三维图
persp(x = 1:g[1], # 代数范围
      y = possible, # 可能的等位基因数量
      z = Prob, # 概率矩阵
      theta = 60, phi = 20, # 视图旋转角度
      xlab = 'Generations', ylab = 'Number of alleles', # x轴和y轴标签
      zlab = 'Probability', # z轴标签
      shade = 0.3 # 表面阴影强度
)

# 加载popgenr包中的果蝇数据集
library('popgenr')
data(fly)

# 使用相同的参数绘制实际观测数据的图
persp(x = 1:g[1],
      y = possible,
      z = fly[-1,], 
      theta = 60, phi = 20,
      xlab = 'Generations', ylab = 'Number of alleles',
      zlab = 'Observations', # z轴标签改为'Observations'
      shade = 0.3
)

观察这两个图，我们应该看到bw75等位基因的行为与我们假设的A等位基因的二项式抽样模型非常相似。在模拟和实验结束时，我们看到最高峰出现在最左侧的“灭绝”状态（零个A/bw75等位基因）和最右侧的“固定”状态（与位点一样多的A/bw75等位基因）。这表明，如果没有其他力量作用于单个等位基因，随着时间的推移，每一个新的变体要么消失，要么不再是明显的变体，因为它成为了位点上的唯一等位基因。中央峰在广泛的概率范围内的衰减应该让人想起物理学中的热扩散（激发了中性理论的“扩散近似”）；然而，频率为零和一的吸收边界最终是概率集中的地方。那么为什么我们应该在种群内观察到任何水平的遗传变异呢？突变率当然是一个因素，但我们已经看到，新的变体几乎在出现时就很有可能消失。所以，让我们考虑一些原因，为什么看似中性的等位基因会在种群中持续存在。我们之前看到种群大小可以影响单个等位基因的固定和灭绝概率。让我们再次可视化这种效应，但这次同时观察多个不同的模拟。我们将根据种群大小和它们的起始频率，使用rbinom随机抽取每一代的等位基因数量：

我们分别对于10个的群体、100个群体及1000个群体进行模拟：

par(mfrow=c(2,2))
# 初始化等位基因频率
init_p <- 0.05

# 设定模拟的代数
gen <- 100

# 设定重复实验的次数
reps <- 10

# 生成一系列颜色，用于区分不同的重复实验
colors <- rainbow(reps)

# 设定群体大小
N <- 100

# 创建一个空白的绘图区域，设置x轴和y轴的范围及标签
plot(x=NULL, y=NULL, xlim=c(1, gen), ylim=c(0,1),
     xlab="Generations", ylab="Allele frequency")

# 对于每次重复实验
for(i in 1:reps){
  # 初始化等位基因频率
  p <- init_p
  # 对于每一代（除了第一代）
  for(j in 1:(gen-1)){
    # 使用二项分布模拟等位基因的传递
    a <- rbinom(n=1, size=2*N, prob=p[j])
    # 计算新的等位基因频率
    f <- a/(2*N)
    # 更新等位基因频率向量
    p <- c(p,f)
  }
  # 在图中画出这一重复实验的等位基因频率变化轨迹
  lines(x=1:gen, y=p, lwd=2, col=colors[i])
}

随着种群规模的增大，等位基因似乎会在更长的时间内持续存在。然而，也尝试改变这个模拟中的世代数（Gen）。生物体的世代时间可能有很大的不同：在我们人类漫长一代的时间里，你可能会看到几百代、几千代或更多的果蝇。如果你让它们运行1,000代，而不是100代，较大的种群与较小的种群相比如何？随着时间的推移，变异的丧失几乎是肯定的，但种群大小和世代时间等因素显著影响我们预期等位基因在种群中丧失的速度。最后，这里有一个遗传漂变的模拟，它跟踪不同的轨迹，并在完成时对它们进行平均。你可以调整起始频率和种群大小，虽然个体轨迹遍布各处，但最终用粗黑线绘制的平均值基本保持不变。

# 初始化等位基因频率
init_p <- 0.25

# 设定模拟的代数
gen <- 100

# 设定重复实验的次数
reps <- 500

# 生成一系列颜色，用于区分不同的重复实验
colors <- rainbow(reps)

# 设定群体大小
N <- 100

# 创建一个空白的绘图区域，设置x轴和y轴的范围及标签
plot(x=NULL, y=NULL, xlim=c(1, gen), ylim=c(0,1),
     xlab="Generations", ylab="Allele frequency")

# 初始化一个空矩阵Freq，用于保存每次重复实验的等位基因频率
Freq <- NULL

# 对于每次重复实验
for(i in 1:reps){
  # 初始化等位基因频率
  p <- init_p
  # 对于每一代（除了第一代）
  for(j in 1:(gen-1)){
    # 使用二项分布模拟等位基因的传递
    a <- rbinom(n=1, size=2*N, prob=p[j])
    # 计算新的等位基因频率
    f <- a/(2*N)
    # 更新等位基因频率向量
    p <- c(p,f)
  }
  # 将这次重复实验的等位基因频率保存到Freq矩阵中
  Freq <- rbind(Freq, p)
  # 在图中画出这一重复实验的等位基因频率变化轨迹
  lines(x=1:gen, y=p, lwd=2, col=colors[i])
}

# 计算所有重复实验的平均等位基因频率，并在图中以黑色线条画出
lines(1:gen, colMeans(Freq), lwd=2, col="black")