#An example of several scripts I used to calculate fd on simulated sequence data with no gene flow. I may be talking to the void, but my apologies in advance if you attempt to use parts of it and it doesn't work out too well for you the first time (or several) around, I am a noob at shareable code. ######################### #Generate 100000 sequences under model of no gene flow but population changes over time using ms (and MSMC estimates of Ne changes: FigS13-S14) REPS=100000 for I in $(seq 1 $REPS) do #faster mutation rate simulation Ne changes #/proj/cmarlab/users/emilie/software/msmove-master/msmove 8 1 -I 4 2 1 3 2 0.000001 -g 1 40 -g 2 40 -g 3 1.59 -g 4 3.9 -eg 0.025 1 0 -eg 0.025 2 0 -eg 0.18 3 0 -eg 0.1 4 0 -en 0.25 1 0.1 -en 0.25 2 0.1 -en 0.6 3 0.5 -en 1.25 4 0.5 -T | tail -n +3 | head -n 1 | grep -v // > ../treefiles/treefile.$I.out #slower mutation rate simulation Ne changes /proj/cmarlab/users/emilie/software/msmove-master/msmove 8 1 -I 4 2 1 3 2 0.000001 -g 1 40 -g 2 40 -g 3 1.59 -g 4 3.9 -eg 0.25 1 0 -eg 0.25 2 0 -eg 1.8 3 0 -eg 1 4 0 -en 2.5 1 0.1 -en 2.5 2 0.1 -en 6 3 0.5 -en 12.5 4 0.5 -T | tail -n +3 | head -n 1 | grep -v // > ../treefiles/treefile.$I.out /proj/cmarlab/users/emilie/software/SeqGen/seq-gen -mGTR -l 2000 < ../treefiles/treefile.$I.out > ../seqfiles/seqfile.$I.out done ######################## ######################## #concatenate and reformat simulated sequences to prep for conversion into fd format sed -i '1d' ./seqfiles/seqfile.1.out sort -k1 -n -o ./seqfiles/seqfile.1.out ./seqfiles/seqfile.1.out #sed -i -e 's/^[0-9]\+\s\+//g' ms.geneflowsim.seqfile.1.out for I in {2..100000}; do sed -i '1d' ./seqfiles/seqfile.$I.out sort -k1 -n -o ./seqfiles/seqfile.$I.out ./seqfiles/seqfile.$I.out sed -i -e 's/^[0-9]\+\s\+//g' ./seqfiles/seqfile.$I.out done paste -d '' ./seqfiles/seqfile.{1..50000}.out > nogeneflow_catseqfiles1.test paste -d '' ./seqfiles/seqfile.{50001..100000}.out > nogeneflow_catseqfiles2.test paste -d '' ./nogeneflow_catseqfiles1.test ./nogeneflow_catseqfiles2.test > cat_seq_no_geneflowsims_allNe_slowerm.out sed -i '1i7 200000000' cat_seq_no_geneflowsims.out ######################## ######################## #Generate the .geno file required to calculate fd using the ABBABABAwindows.py script and then run fd in windows on simulated sequences python /proj/cmarlab/users/emilie/software/genomics_general-master/seqToGeno.py -s ../cat_seq_no_geneflowsims_allNe_slowerm.out -f phylip -g cat_seq_no_geneflowsims_allNe_slowerm.geno python /proj/cmarlab/users/emilie/software/genomics_general-master/ABBABABAwindows.py -g /pine/scr/e/j/ejr/Cichlids/seq_sims/scripts/cat_seq_no_geneflowsims_allNe_slowerm.geno -o /pine/scr/e/j/ejr/Cichlids/data_files/cichlid_fd_no_geneflowsim_2k_allNe_slowerm.csv -f diplo --windType coordinate -w 2000 -P1 P1 1,2 -P2 P2 3 -P3 P3 4,5,6 -O Out 7,8