1

我正在尝试在 R 中运行wright-fisher遗传漂移模型的模拟。

# 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 轴的长度）... 在此处输入图像描述

4

2 回答 2

3

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

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

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)

在此处输入图像描述

于 2013-10-10T12:06:02.250 回答
0 
# 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()
于 2016-03-31T16:39:14.677 回答