Post date: Dec 09, 2015 8:2:55 PM
I am testing popanc on simulations with migration as requested by a reviewer. All files are in /labs/evolution/projects/popanc_sims/sims/ or subdirectories of this directory.
I first used admix to simulate neutral admixture with on-going gene flow. Simulations were for 200 generations with m = 0.005 or m = 0.05 (these numbers give the total proportion of migrants). Other than that, conditions match those described here. The perl script used to run admix is called runAdmixPopAnc3.pl and is pasted below:
$nmarkersChr = 10001; ## number of markers per LG
makenames();
$model = 'underdom'; ## form of selection: underdom or epistasis
$multadd = 'mult'; ## fitness with multiple loci is multiplicative 'mult' or additive
$gen = 0; ## no. of generations
$rep = 10; ## no. of replicates
$pop = 500; ## population size = sample size?
$mig = 0.0; ## migration rate
$nsel= 1; ## total no. of selected loci
$nselperchr= 1; ## no. of selected loci on each chromosome
$distselloci = 1; ## spacing between selected loci, if more than one per chrom.
$sel = 0.0; ## selection coefficient
$gen = 200; ##
foreach $mig (0.005, 0.05){
$outfile = "admix_demog_mig".$mig.".txt";
print "$outfile: $gen $sel $mig $pop\n";
open(OUTFILE, "> $outfile") or die "Couldn't open outfile\n";
print OUTFILE "model,stage,multadd,s,nselloci,nsellociperchr,distselloci,mig,gen,rep,ind,h,het,$loci\n";
close(OUTFILE);
system "admix $model viability $sel $nsel $nselperchr $distselloci $multadd $gen $rep $pop $mig $nmarkersChr $outfile";
}
sub makenames {
foreach $j (1..10){
foreach $i (1..$nmarkersChr){
push @loci, "c$j.$i";
}
}
$loci = join ",", @loci;
}
Next, I used simGenotypesFstMig.pl to generate genotypes from the ancestry data (as previously described), except that here I explicitly modelled drift with gene flow. I than used reHeader.pl to add headers to these files.
Finally, I ran popanc on the data sets. As with other simulated data sets, I ran 30,000 iterations with a 10,000 iteration burn-in and thinning interval of 5. I ran two chains per analysis and used windows of 4, 10, or 20 SNPs. Results will be in /labs/evolution/projects/popanc_sims/mcmc/. Here is an example command:
cd /local/scratch/
sleep 5
popanc -o estpanc_n20migMig0.05Rep9Chain1.hdf5 -m 30000 -b 10000 -t 5 -n 20 -d 20 -s 1 -u 10 -l 0.1 -a 0.1 -z 1 /labs/evolution/projects/popanc_sims/sims/mig/genoP0F_r9gen_demog_mig0.05.txt /labs/evolution/projects/popanc_sims/sims/mig/genoP1F_r9gen_demog_mig0.05.txt /labs/evolution/projects/popanc_sims/sims/mig/genoAdmxF_r9gen_demog_mig0.05.txt
scp estpanc_n20migMig0.05Rep9Chain1.hdf5 /labs/evolution/projects/popanc_sims/mcmc/