#admixture simulation with small Ne, explosive breeding, and nonoverlapping generations library(stats) library(genetics) parents=genotype(a1=c(rep("A",5),rep("a",5)),a2=c(rep("A",5),rep("a",5))) kids.fn=function(parents,F){ if(prod(parents)==1) kids=c("AA"=F,"Aa"=0,"aa"=0) if(prod(parents)==2) kids=as.vector(rmultinom(1,F,p=c(.5,.5,0))) if(prod(parents)==3) kids=c("AA"=0,"Aa"=F,"aa"=0) if(prod(parents)==4) kids=as.vector(rmultinom(1,F,p=c(.25,.5,.25))) if(prod(parents)==6) kids=as.vector(rmultinom(1,F,p=c(0,.5,.5))) if(prod(parents)==9) kids=c("AA"=0,"Aa"=0,"aa"=F) kids } Np=5 F=1000 n=41 L=64 T=80 genotypes=genotype(c("A/A","A/a","a/a")) Freqs=matrix(NA,ncol=L,nrow=T) Fis=matrix(NA,ncol=L,nrow=T) for(j in 1:L){ Geno=c("AA"=1200,"Aa"=0,"aa"=800) Geno=rbind(Geno,Geno,Geno) for(i in 3:T){ probs=colSums(Geno[c(i-2,i-1,i),]) pairs=cbind(sample(genotypes,size=Np,replace=TRUE,prob=probs),sample(genotypes,size=Np,replace=TRUE,prob=probs)) kids=c("AA"=0,"Aa"=0,"aa"=0) for(m in 1:Np){ kids=kids+kids.fn(pairs[m,],F) } Geno=rbind(Geno,kids) } for(i in 1:T){ Sam=sample(genotypes,size=n,replace=TRUE,prob=Geno[i,]) Freqs[i,j]=summary(Sam)\$allele.freq[1,1] Fis[i,j]=diseq(Sam,R=0)\$r.overall } } Pbar=rowMeans(Freqs/(2*n)) Pvar=apply(Freqs/(2*n),1,var) Fls=Pvar/(Pbar*(1-Pbar)) Pfix=rowSums(Freqs==0)+rowSums(Freqs==2*n) Pi=(Freqs+1)/(2*n+2) Vmi=matrix(NA,ncol=L,nrow=T) for(i in 1:T){ for(j in 1:L){ Vmi[i,j]=Pi[i,j]*(1-Pi[i,j])/(1/(2*n)+Fls[i]*(1-1/(2*n))) }} LongX2=rowSums((Pi-Pbar)^2/Vmi) cbind(Fls,Pfix,LongX2,P=pchisq(LongX2,df=L-1,lower.tail=FALSE)) ####### ######## ######## For a range of Np N=c(2,4,8,16,32,64) MeF=Mefix=MeLong=matrix(NA,ncol=6,nrow=44) for(z in 1:6){ Np=N[z] F=1000 n=55 L=64 T=80 genotypes=genotype(c("A/A","A/a","a/a")) Freqs=matrix(NA,ncol=L,nrow=T) Fis=matrix(NA,ncol=L,nrow=T) for(j in 1:L){ Geno=c("AA"=246,"Aa"=0,"aa"=1754) Geno=rbind(Geno,Geno,Geno) for(i in 3:T){ probs=colSums(Geno[c(i-2,i-1,i),]) pairs=cbind(sample(genotypes,size=Np,replace=TRUE,prob=probs),sample(genotypes,size=Np,replace=TRUE,prob=probs)) kids=c("AA"=0,"Aa"=0,"aa"=0) for(m in 1:Np){ kids=kids+kids.fn(pairs[m,],F) } Geno=rbind(Geno,kids) } for(i in 1:T){ Sam=sample(genotypes,size=n,replace=TRUE,prob=Geno[i,]) Freqs[i,j]=summary(Sam)\$allele.freq["A",1] Fis[i,j]=diseq(Sam,R=0)\$r.overall } } Pbar=rowMeans(Freqs/(2*n)) Pvar=apply(Freqs/(2*n),1,var) Fls=Pvar/(Pbar*(1-Pbar)) Pfix=rowSums(Freqs==2*n) Pi=(Freqs+1)/(2*n+2) Vmi=matrix(NA,ncol=L,nrow=T) for(i in 1:T){ for(j in 1:L){ Vmi[i,j]=Pi[i,j]*(1-Pi[i,j])/(1/(2*n)+Fls[i]*(1-1/(2*n))) }} LongX2=rowSums((Pi-Pbar)^2/Vmi) MeLong[,z]=LongX2 MeF[,z]=Fls Mefix[,z]=Pfix } quartz(height=12) par(mfrow=c(3,1)) plot(MeF,MeLong) points(.548,854.5,pch=4) plot(Mefix,MeLong) points(3,854.5,pch=4) plot(MeF,Mefix) points(.548,3,pch=4) ###### Long.sim=function(Np,F,n,L,T,reps=5,M){ FLS=fix=Long=matrix(NA,ncol=reps,nrow=T) X=round(M*2000) Y=round((1-M)*2000) for(z in 1:reps){ genotypes=genotype(c("A/A","A/a","a/a")) Freqs=matrix(NA,ncol=L,nrow=T) Fis=matrix(NA,ncol=L,nrow=T) for(j in 1:L){ Geno=c("AA"=X,"Aa"=0,"aa"=Y) Geno=rbind(Geno,Geno,Geno) for(i in 3:T){ probs=colSums(Geno[c(i-2,i-1,i),]) pairs=cbind(sample(genotypes,size=Np,replace=TRUE,prob=probs),sample(genotypes,size=Np,replace=TRUE,prob=probs)) kids=c("AA"=0,"Aa"=0,"aa"=0) for(m in 1:Np){ kids=kids+kids.fn(pairs[m,],F=F) } Geno=rbind(Geno,kids) } for(i in 1:T){ Sam=sample(genotypes,size=n,replace=TRUE,prob=Geno[i,]) Freqs[i,j]=summary(Sam)\$allele.freq["A",1] Fis[i,j]=diseq(Sam,R=0)\$r.overall } } Pbar=rowMeans(Freqs/(2*n)) Pvar=apply(Freqs/(2*n),1,var) Fls=Pvar/(Pbar*(1-Pbar)) Pfix=rowSums(Freqs==2*n) Pi=(Freqs+1)/(2*n+2) Vmi=matrix(NA,ncol=L,nrow=T) for(i in 1:T){ for(j in 1:L){ Vmi[i,j]=Pi[i,j]*(1-Pi[i,j])/(1/(2*n)+Fls[i]*(1-1/(2*n))) }} LongX2=rowSums((Pi-Pbar)^2/Vmi) Long[,z]=LongX2 FLS[,z]=Fls fix[,z]=Pfix } list(FLS,fix,Long) } BS2=Long.sim(Np=2,F=1000,n=41,L=64,T=80,reps=500,M=.607) quartz(width=10,height=10) par(mfrow=c(5,5)) plot(BS2[[1]],BS2[[3]],xlab="FLS",ylab="X2",main="BS2",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.073,543.6,pch=4,cex=2) BS4=Long.sim(Np=4,F=1000,n=41,L=64,T=80,reps=500,M=.607) plot(BS4[[1]],BS4[[3]],xlab="FLS",ylab="X2",main="BS4",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.073,543.6,pch=4,cex=2) BS8=Long.sim(Np=8,F=1000,n=41,L=64,T=80,reps=500,M=.607) plot(BS8[[1]],BS8[[3]],xlab="FLS",ylab="X2",main="BS8",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.073,543.6,pch=4,cex=2) BS16=Long.sim(Np=16,F=1000,n=41,L=64,T=80,reps=500,M=.607) plot(BS16[[1]],BS16[[3]],xlab="FLS",ylab="X2",main="BS16",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.073,543.6,pch=4,cex=2) BS32=Long.sim(Np=32,F=1000,n=41,L=64,T=80,reps=500,M=.607) plot(BS32[[1]],BS32[[3]],xlab="FLS",ylab="X2",main="BS32",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.073,543.6,pch=4,cex=2) ME2=Long.sim(Np=2,F=1000,n=55,L=64,T=80,reps=500,M=.123) plot(ME2[[1]],ME2[[3]],xlab="FLS",ylab="X2",main="ME2",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.548,854.5,pch=4,cex=2) ME4=Long.sim(Np=4,F=1000,n=55,L=64,T=80,reps=500,M=.123) plot(ME4[[1]],ME4[[3]],xlab="FLS",ylab="X2",main="ME4",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.548,854.5,pch=4,cex=2) ME8=Long.sim(Np=8,F=1000,n=55,L=64,T=80,reps=500,M=.123) plot(ME8[[1]],ME8[[3]],xlab="FLS",ylab="X2",main="ME8",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.548,854.5,pch=4,cex=2) ME16=Long.sim(Np=16,F=1000,n=55,L=64,T=80,reps=500,M=.123) plot(ME16[[1]],ME16[[3]],xlab="FLS",ylab="X2",main="ME16",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.548,854.5,pch=4,cex=2) ME32=Long.sim(Np=32,F=1000,n=55,L=64,T=80,reps=500,M=.123) plot(ME32[[1]],ME32[[3]],xlab="FLS",ylab="X2",main="ME32",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.548,854.5,pch=4,cex=2) H2=Long.sim(Np=2,F=1000,n=51,L=64,T=80,reps=500,M=.642) plot(H2[[1]],H2[[3]],xlab="FLS",ylab="X2",main="H2",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.142,663.3,pch=4,cex=2) H4=Long.sim(Np=4,F=1000,n=51,L=64,T=80,reps=500,M=.642) plot(H4[[1]],H4[[3]],xlab="FLS",ylab="X2",main="H4",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.142,663.3,pch=4,cex=2) H8=Long.sim(Np=8,F=1000,n=51,L=64,T=80,reps=500,M=.642) plot(H8[[1]],H8[[3]],xlab="FLS",ylab="X2",main="H8",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.142,663.3,pch=4,cex=2) H16=Long.sim(Np=16,F=1000,n=51,L=64,T=80,reps=500,M=.642) plot(H16[[1]],H16[[3]],xlab="FLS",ylab="X2",main="H16",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.142,663.3,pch=4,cex=2) H32=Long.sim(Np=32,F=1000,n=51,L=64,T=80,reps=500,M=.642) plot(H32[[1]],H32[[3]],xlab="FLS",ylab="X2",main="H32",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.142,663.3,pch=4,cex=2) S2=Long.sim(Np=2,F=1000,n=56,L=64,T=80,reps=500,M=.340) plot(S2[[1]],S2[[3]],xlab="FLS",ylab="X2",main="S2",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.161,1088.7,pch=4,cex=2) S4=Long.sim(Np=4,F=1000,n=56,L=64,T=80,reps=500,M=.340) plot(S4[[1]],S4[[3]],xlab="FLS",ylab="X2",main="S4",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.161,1088.7,pch=4,cex=2) S8=Long.sim(Np=8,F=1000,n=56,L=64,T=80,reps=500,M=.340) plot(S8[[1]],S8[[3]],xlab="FLS",ylab="X2",main="S8",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.161,1088.7,pch=4,cex=2) S16=Long.sim(Np=16,F=1000,n=56,L=64,T=80,reps=500,M=.340) plot(S16[[1]],S16[[3]],xlab="FLS",ylab="X2",main="S16",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.161,1088.7,pch=4,cex=2) S32=Long.sim(Np=32,F=1000,n=56,L=64,T=80,reps=500,M=.340) plot(S32[[1]],S32[[3]],xlab="FLS",ylab="X2",main="S32",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.161,1088.7,pch=4,cex=2) T2=Long.sim(Np=2,F=1000,n=52,L=64,T=80,reps=500,M=.545) plot(T2[[1]],T2[[3]],xlab="FLS",ylab="X2",main="T2",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.181,755.2,pch=4,cex=2) T4=Long.sim(Np=4,F=1000,n=52,L=64,T=80,reps=500,M=.545) plot(T4[[1]],T4[[3]],xlab="FLS",ylab="X2",main="T4",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.181,755.2,pch=4,cex=2) T8=Long.sim(Np=8,F=1000,n=52,L=64,T=80,reps=500,M=.545) plot(T8[[1]],T8[[3]],xlab="FLS",ylab="X2",main="T8",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.181,755.2,pch=4,cex=2) T16=Long.sim(Np=16,F=1000,n=52,L=64,T=80,reps=500,M=.545) plot(T16[[1]],T16[[3]],xlab="FLS",ylab="X2",main="T16",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.181,755.2,pch=4,cex=2) T32=Long.sim(Np=32,F=1000,n=52,L=64,T=80,reps=500,M=.545) plot(T32[[1]],T32[[3]],xlab="FLS",ylab="X2",main="T32",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.181,755.2,pch=4,cex=2) #write T8=read("T8.txt") T16=read.table("T16.txt") T32=read.table("T32.txt") T8=list(as.matrix(T8[,1:500]),as.matrix(T8[,501:1000]),as.matrix(T8[,1001:1500])) T16=list(as.matrix(T16[,1:500]),as.matrix(T16[,501:1000]),as.matrix(T16[,1001:1500])) T32=list(as.matrix(T32[,1:500]),as.matrix(T32[,501:1000]),as.matrix(T32[,1001:1500])) quartz(width=10,height=10) par(mfrow=c(5,5)) plot(BS2[[1]],BS2[[3]],xlab="FLS",ylab="X2",main="BS2",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.073,543.6,pch=4,cex=2) plot(BS4[[1]],BS4[[3]],xlab="FLS",ylab="X2",main="BS4",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.073,543.6,pch=4,cex=2) plot(BS8[[1]],BS8[[3]],xlab="FLS",ylab="X2",main="BS8",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.073,543.6,pch=4,cex=2) plot(BS16[[1]],BS16[[3]],xlab="FLS",ylab="X2",main="BS16",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.073,543.6,pch=4,cex=2) plot(BS32[[1]],BS32[[3]],xlab="FLS",ylab="X2",main="BS32",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.073,543.6,pch=4,cex=2) plot(ME2[[1]],ME2[[3]],xlab="FLS",ylab="X2",main="ME2",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.548,854.5,pch=4,cex=2) plot(ME4[[1]],ME4[[3]],xlab="FLS",ylab="X2",main="ME4",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.548,854.5,pch=4,cex=2) plot(ME8[[1]],ME8[[3]],xlab="FLS",ylab="X2",main="ME8",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.548,854.5,pch=4,cex=2) plot(ME16[[1]],ME16[[3]],xlab="FLS",ylab="X2",main="ME16",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.548,854.5,pch=4,cex=2) plot(ME32[[1]],ME32[[3]],xlab="FLS",ylab="X2",main="ME32",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.548,854.5,pch=4,cex=2) plot(H2[[1]],H2[[3]],xlab="FLS",ylab="X2",main="H2",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.142,663.3,pch=4,cex=2) plot(H4[[1]],H4[[3]],xlab="FLS",ylab="X2",main="H4",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.142,663.3,pch=4,cex=2) plot(H8[[1]],H8[[3]],xlab="FLS",ylab="X2",main="H8",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.142,663.3,pch=4,cex=2) plot(H16[[1]],H16[[3]],xlab="FLS",ylab="X2",main="H16",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.142,663.3,pch=4,cex=2) plot(H32[[1]],H32[[3]],xlab="FLS",ylab="X2",main="H32",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.142,663.3,pch=4,cex=2) plot(S2[[1]],S2[[3]],xlab="FLS",ylab="X2",main="S2",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.161,1088.7,pch=4,cex=2) plot(S4[[1]],S4[[3]],xlab="FLS",ylab="X2",main="S4",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.161,1088.7,pch=4,cex=2) plot(S8[[1]],S8[[3]],xlab="FLS",ylab="X2",main="S8",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.161,1088.7,pch=4,cex=2) plot(S16[[1]],S16[[3]],xlab="FLS",ylab="X2",main="S16",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.161,1088.7,pch=4,cex=2) plot(S32[[1]],S32[[3]],xlab="FLS",ylab="X2",main="S32",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.161,1088.7,pch=4,cex=2) plot(T2[[1]],T2[[3]],xlab="FLS",ylab="X2",main="T2",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.181,755.2,pch=4,cex=2) plot(T4[[1]],T4[[3]],xlab="FLS",ylab="X2",main="T4",ylim=c(0,1200),xlim=c(0,1),pch=16,cex=.3) points(.181,755.2,pch=4,cex=2) plot(T8[[1]],T8[[3]],xlab="FLS",ylab="X2",main="T8",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.181,755.2,pch=4,cex=2) plot(T16[[1]],T16[[3]],xlab="FLS",ylab="X2",main="T16",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.181,755.2,pch=4,cex=2) plot(T32[[1]],T32[[3]],xlab="FLS",ylab="X2",main="T32",ylim=c(0, 1200),xlim=c(0,1),pch=16,cex=.3) points(.181,755.2,pch=4,cex=2) ### FLS over time quartz(width=10,height=12) par(mfrow=c(5,5)) plot(rep(1:80,500),BS2[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="BS2") abline(h=0.073,lty=2) plot(rep(1:80,500),BS4[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="BS4") abline(h=0.073,lty=2) plot(rep(1:80,500),BS8[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="BS8") abline(h=0.073,lty=2) plot(rep(1:80,500),BS16[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="BS16") abline(h=0.073,lty=2) plot(rep(1:80,500),BS32[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="BS32") abline(h=0.073,lty=2) plot(rep(1:80,500),ME2[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="ME2") abline(h=0.548,lty=2) plot(rep(1:80,500),ME4[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="ME8") abline(h=0.548,lty=2) plot(rep(1:80,500),ME8[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="ME4") abline(h=0.548,lty=2) plot(rep(1:80,500),ME16[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="ME16") abline(h=0.548,lty=2) plot(rep(1:80,500),ME32[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="ME32") abline(h=0.548,lty=2) plot(rep(1:80,500),H2[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="H2") abline(h=0.142,lty=2) plot(rep(1:80,500),H4[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="H4") abline(h=0.142,lty=2) plot(rep(1:80,500),H8[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="H8") abline(h=0.142,lty=2) plot(rep(1:80,500),H16[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="H16") abline(h=0.142,lty=2) plot(rep(1:80,500),H32[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="H32") abline(h=0.142,lty=2) plot(rep(1:80,500),S2[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="S2") abline(h=0.161,lty=2) plot(rep(1:80,500),S4[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="S4") abline(h=0.161,lty=2) plot(rep(1:80,500),S8[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="S8") abline(h=0.161,lty=2) plot(rep(1:80,500),S16[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="S16") abline(h=0.161,lty=2) plot(rep(1:80,500),S32[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="S32") abline(h=0.161,lty=2) plot(rep(1:80,500),T2[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="T2") abline(h=0.181,lty=2) plot(rep(1:80,500),T4[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="T4") abline(h=0.181,lty=2) plot(rep(1:80,500),T8[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="T8") abline(h=0.181,lty=2) plot(rep(1:80,500),T16[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="T16") abline(h=0.181,lty=2) plot(rep(1:80,500),T32[[1]],xlab="Years",ylab="FLS",ylim=c(0,1),main="T32") abline(h=0.181,lty=2)