Migration
Validation of migration
MimicrEE2 supports migration (in w-mode, qt-mode and qff-mode).
To validate the correct implementation of migration in MimicrEE2 we run 100 independent simulations of a scenario where we have a population of size N=1000 individuals, that is evolving for 60 generations. No beneficial SNPs are specified therefore neutral simulations are performed. The initial allele frequency of the observed SNP is 0.1
Every generation 10% of the individuals from a source population are migrating to the evolved population, replacing 10% of its individuals. In the source population the allele frequency of the observed SNP is 1 (the allele is fixed).
We thus expect that the frequency of the observed SNP increases in the evolving population unit it eventually becomes fixed. We compared the observed trajectory of the observed SNPs with the expected trajectory (see Material and Methods for the theoretical expectation).
The grey lines correspond to simulations with MimicrEE2, the dark line to the theoretical expectation and the dashed line to the allele frequency in the source population. As expected the allele frequency in the receiver population increases until the allele becomes fixed. This shows that the theoretical and observed allele frequency changes due to migration are matching.
Materials and methods
All input files and scripts used for the simulations can be found here.
For more details on the input files see input files.
First, 100 independent simulations were performed with MimicrEE2. Output is requested all 5 generations.
:::bash
for i in {1..100}; do java -jar ../mim2.jar qt --haplotypes-g0 basepop --recombination-rate one_locus_rr --effect-size effect_size --ve 0.0 --output-sync test_generations_${i}.txt --selection-regime selection_regime --threads 6 --snapshots 0,5,10,15,20,25,30,35,40,45,50,55,60 --migration-regime migration_rate; done
After the simulations are finished we can convert the sync files to allele frequency files with the following python script:
:::bash
for i in {1..100}; do python freq_validation_more_than_one.py --n 1 --input test_generations_${i}.txt --output test_generations_${i}.fq --chosen-snp 100 --minor-allele 0 --generations-number 12; done
##concatenate output files
cat test_generations_{1..100}.fq > freq_migration
Finally we run a function to plot the simulated and expected allele frequences:
:::R
#p0 is the initial allele frequency of the receiver population
#pm0 is the allele frequency in the migrant population
#t is the number of generations
# m is the migration rate
#repl is the number of replicate runs form mimicrEE2
migration<- function(p0,pm0, t, m, repl ){
##set your current working directory
setwd("my/working/dir")
traj<- matrix(NA, ncol=3, nrow=max(t)+1)
trajm<- matrix(NA, ncol=3, nrow=max(t)+1)
traj[1,1]<-p0
trajm[,1]<-pm0
traj[1,2]<-0
trajm[1,2]<-0
g<-1
p<-p0
pm<-pm0
while(g <= max(t)){
p<-(1-m)*p + m*pm
traj[g+1,1]<-p
traj[g+1,2]<-g
trajm[g+1,2]<-g
g<-g+1
}
traj<- as.data.frame(traj)
trajm<- as.data.frame(trajm)
traj[,3]<- "black"
trajm[,3]<- "black"
colnames(traj)<- c("freq", "generations", "color")
colnames(trajm)<- c("freq", "generations", "color")
print(traj)
print(trajm)
generations=t
mimicree_migration<- read.table("freq_migration")
##5 means that we get estimates after every five generations (e.g for generation 0,5,10,15..)
mimicree_migration[,2]<- rep(seq(0,generations,5), times=repl)
##13 is the count of generations for which estimates are provided
##(e.g if we have values for generations 0,5,10,15,20,25, then this value should be equal to 6)
mimicree_migration[,3]<- rep(c(1:repl),each=13)
mimicree_migration[,3]<- as.character(mimicree_migration[,3])
colnames(mimicree_migration)<- colnames(traj)
library(ggplot2)
ggplot(rbind(mimicree_migration, traj), aes(x=generations, y=freq, color=color,size=color)) + geom_line()+
geom_line(data=trajm,aes(x=generations, y=freq), color="black",size=0.5, linetype="dashed")+
theme_bw() +
labs(x="Generations", y="Allele Frequency(%)") +
scale_color_manual(values=c(rep("grey",times=repl),"black"),guide=FALSE) +
scale_size_manual(values=c(rep(0.5,times=repl),0.7),guide=FALSE)
}
##run function
migration(0.1, 1, 60, 0.1, 100)
