Post date: Mar 22, 2015 4:55:8 PM
I modified simGenotypesFst.pl to extract the ancestry information and prepare genotype files for popanc. I made a few changes. Most notably, I only grab 20,002 SNPs (first two chromosomes), allele frequencies for each source population are U(0.05,0.95) with the two populations being independent, and I automatically add one heterozygote to P0 if both populations are fixed for the same allele (this should be very very rare, but it is important to avoid having a singular matrix). After this the reHeader.pl script is used to add the header. Then its time to let the program rip.
Here is a version of the R script it writes for one replicate:
g<-G[G[,1]==10,-c(1:4)]
gAdmx<-t(g[1:50,])
gP0<-matrix(NA,nrow=20002,ncol=50)
gP1<-gP0
pP0<-runif(20002,0.05,0.95)
pP1<-runif(20002,0.05,0.95)
for(i in 1:20002){
for(j in 1:50){
gP0[i,j]<-rbinom(1,2,pP0[i])
gP1[i,j]<-rbinom(1,2,pP1[i])
if((sum(gP0[i,]) + sum(gP1[i,])) == 0){gP0[i,1]<-1}
if((sum(gP0[i,]) + sum(gP1[i,])) == 2){gP0[i,1]<-1}
}
}
for(i in 1:20002){
t<-20/500
f<-exp(-t) * (exp(t) - 1)
S<--1 * (f-1)/f
if(pP0[i] > 0 & pP0[i] < 1){pP0[i]<-rbeta(1,pP0[i] * S,(1 - pP0[i]) * S)}
if(pP1[i] > 0 & pP1[i] < 1){pP1[i]<-rbeta(1,pP1[i] * S,(1 - pP1[i]) * S)}
for(j in 1:50){
if(gAdmx[i,j] == 0) {gAdmx[i,j]<-rbinom(1,2,pP0[i])}
else if(gAdmx[i,j] == 2) {gAdmx[i,j]<-rbinom(1,2,pP1[i])}
else if(gAdmx[i,j] == 1) {gAdmx[i,j]<-rbinom(1,1,pP0[i]) + rbinom(1,1,pP1[i])}
}
}
write.table(cbind(llist,gAdmx),file="genoAdmxF_r10gen_demog_gens20.txt",row.names=F,col.names=F,quote=F)
write.table(cbind(llist,gP0),file="genoP0F_r10gen_demog_gens20.txt",row.names=F,col.names=F,quote=F,)
write.table(cbind(llist,gP1),file="genoP1F_r10gen_demog_gens20.txt",row.names=F,col.names=F,quote=F,)