使用 R 进行遗传漂移的 Wright-Fisher 模拟答案

【问题标题】：Wright-Fisher Simulation of Genetic Drift using R使用 R 进行遗传漂移的 Wright-Fisher 模拟
【发布时间】：2013-10-10 11:25:19
【问题描述】：

我正在尝试在 R 中模拟wright-fisher model 的遗传漂变。

# Wright-Fisher simulation
# n = number of individuals
# f = number of focal alleles at base population
n=10
f=1
pop = as.matrix( c( rep(0,n-f), rep(1,f) ) )
pop = as.matrix( sample(pop, n, replace=T) )

这行得通，实际上这是一个复制品，每次我运行最后一行脚本时都是新的一代。我想做但不能做的是有一个循环自动循环 X 代并重复 Y 次重复。

它应该将每一代的结果存储在一个数据框中，然后允许我将它们绘制在一个图表中，看起来像这样（其中 f/n 是等位基因频率，每个复制都由一行，代数决定X轴的长度）...

【问题讨论】：

投票结果很接近，因为“这个问题似乎与帮助中心定义的范围内的编程无关。”你能解释一下你为什么这么认为（这样我可以调整问题以更好地适应）？就我而言，这是一个特定的编程问题。
这可能是因为您对所需的输出有一个想法，但没有充分展示 SO 问题通常需要的必要工作。
假设这是错误的，我尝试了很多东西，昨天花了很多时间阅读但无法弄清楚。

标签： r simulation

【解决方案1】：

这是我几年前写的一个函数。您可以设置弹出大小、要模拟的代数和复制数。

由于您没有显示任何自己的代码，我将由您自己决定如何存储输出。无论如何，这应该会让你继续前进：

Drift_graph = function(t,R){
  N<-250
  p<-0.5
  freq<-as.numeric();
    for( i in 1:t ){
      A1=rbinom(1,2*N,p)
      p=A1/(N*2);
      freq[length(freq)+1]<-p;
    }
  plot(freq,type="l",ylim=c(0,1),col=3,xlab="t",ylab=expression(p(A[1])))
    for(u in 1:R){
      freq1<-as.numeric();
      p<-0.5
        for( j in 1:t ){
          A1=rbinom(1,2*N,p)
          p=A1/(N*2);
          freq1[length(freq1)+1]<-p;
        }
      random<-sample(1:1000,1,replace=F)
      randomcolor<-colors()[random] 
      lines(freq1,type="l",col=(randomcolor))
    }
}

Drift_graph(2000,50)

【讨论】：

感谢您的回答，但它没有将结果输入数据框，这就是我没有点击接受的原因，但我已经投了赞成票。我最终在网上找到了一个模拟器，并根据我的需要对其进行了修改（和注释）。谢谢。

【解决方案2】：

# Pop   = Replicate populations
# Gen   = Generations
# NM    = Male population size
# NF    = Female population size
# P     = Frequency of focal allele

GenDriftSim = function(Pop = Pop, Gen = Gen, NM, NF, P, graph = "y", histo = "y"){
            P = (2*(NM+NF))*P
            NE = round((4*NM*NF)/(NM+NF),0)
            SR = round(NM/NF,2)
            Na = NM+NF
            if(graph=="y"){
            plot(c(0,0),type = "n", main = bquote('N'[M]~'/ N'[F]~'='~.(SR)*', N'[A]~'='~.(Na)*', N'[E]~'='~.(NE)), cex.main = 1, 
                xlim = c(1,Gen), ylim=c(0,1), xlab = "Generations", ylab = "Fequency of focal allele")
            }else{}
            for (i in 1:Pop){
            N = NM+NF
            startA = as.vector(c(rep(1, times = P),rep(0, times = (2*N)-P)))
            Population = matrix(c(
                c(sample(startA, size = 2*N, replace = FALSE)),
                c(rep("M", times = NM), rep("F", times = NF))),
                ncol = 3)
            SimResults[(Gen*i)+1-Gen, 3] <<- sum(as.numeric(Population[,1:2]))/(N*2)
            for(j in 1:(Gen-1)){

                    Population = matrix(c(
                        c(sample(sample(Population[(1:NM),1:2], replace = TRUE),N, replace = TRUE)),
                        c(sample(sample(Population[(1+NM):N,1:2], replace = TRUE),N, replace = TRUE)),
                        c(rep("M", times = NM), rep("F", times = NF))), ncol = 3)
                    SimResults[(Gen*i)+1+j-Gen, 3] <<- sum(as.numeric(Population[,1:2]))/(N*2)
                    }
                s = (i*Gen)-Gen+1; e = i*Gen
                r = as.vector(SimResults[s:e, 3])
                if(graph=="y"){
                points(r~c(1:Gen), type = "l")
                }else{}
            }
            if(histo == "y"){SimResults[,1] = rep(1:Pop, each = Gen)
            SimResults[,2] = rep(1:Gen, times = Pop)
            hist(SimResults[,3][SimResults[,2]==Gen], breaks = 100, cex.lab = 0.7, cex.axis = 0.7, xlim = c(0,1), cex.main = 1, main = bquote('N'[M]~'/ N'[F]~'='~.(SR)*', N'[A]~'='~.(Na)*', N'[E]~'='~.(NE)), xlab = paste0("Frequency of focal allele after ",Gen," Generations"))
            }else{}
}

Pop = 10
Gen = 25
P = 0.5

SimResults = matrix(data = NA, ncol = 3, nrow = Gen*Pop)
GenDriftSim(Pop = Pop, Gen = Gen, NM = 100, NF = 900, P = P, graph = "y",  histo = "n")
GenDriftSim(Pop = Pop, Gen = Gen, NM = 180, NF = 180, P = P, graph = "y",  histo = "n")
dev.off()

【讨论】：