## simple simulation of data = initial, allele freqs, environment, ne, betas ## consider 1000 possible causal variants, but pip = 0.1 for all so expect 100 L<-1000 ## 100 snps N<-10 ## 10 pops. G<-10 ## 10 gens ## sample p from beta, assumes no population structure p<-rbeta(L,.5,.5) ## convert to maf p[p>.5]<-1-p[p>.5] ## write results pmat<-matrix(rep(round(p,4),N),nrow=L,ncol=N,byrow=FALSE) write.table(pmat,file="sim1k_p0.txt",row.names=FALSE,col.names=FALSE,quote=FALSE) ### DO NOT RERUN USE THE SAME AS BEFORE ### ## write Ne, 100 samples from posterior, same beta for each u<-600 l<-400 ne<-matrix(round(rbeta(100*N,10,10) * (u-l) + l,2),nrow=100,ncol=N) write.table(ne,file="sim1k_ne.txt",row.names=FALSE,col.names=FALSE,quote=FALSE) ## random environment from standard normal env<-matrix(round(rnorm((G-1) * N,0,1),3),nrow=G-1,ncol=N) write.table(env,file="sim1k_env.txt",row.names=FALSE,col.names=FALSE,quote=FALSE) ## everything above is held constant, focus on different genetic architectures ## will set vz to 1, so vBV < 1 and is = h2 ## h2 = .3 or .7 ## pip = .1 ## assume no uncertainty in SNPs in model pip<-0.1 sd<-0.12 beta<-rnorm(L,0,sd) vBV<-2 * sum(pip * beta^2 * p * (1-p)) vBV # 0.347069 write.table(cbind(rep(0.1,L),beta),file="sim1k_trait_1.txt",row.names=FALSE,col.names=FALSE,quote=FALSE) sd<-0.17 beta<-rnorm(L,0,sd) vBV<-2 * sum(pip * beta^2 * p * (1-p)) vBV # 0.7066325 write.table(cbind(rep(0.1,L),beta),file="sim1k_trait_2.txt",row.names=FALSE,col.names=FALSE,quote=FALSE)