### written by Kevin Schneider ### ### last edited: 28.07.2020 ### ### Kevin.Schneider@glasgow.ac.uk ### ### KevinSchneider@gmx.at ### require(reshape2) require(ggplot2) require(rcartocolor) require(plyr) require(stringr) require(colorspace) require(gplots) require(pheatmap) require(RColorBrewer) setwd("/Volumes/LaCie/demographic_simulation_paper/fvec_ms_and_vcf_simulation_output_20200709/") path <- "/Volumes/LaCie/demographic_simulation_paper/ms_and_vcf_simulation_output_20200709/" path2 <- "/Volumes/LaCie/demographic_simulation_paper/ms_and_vcf_simulation_output_20200709/rds_dir/" path3 <- "/Users/kevinschneider/Documents/demographic_simulation_paper/analyses_sim20191125_ms_and_vcf_80K/" path4 <- "/Users/kevinschneider/Documents/demographic_simulation_paper/beta_reanalysis/rds_for_Kevin_manhattan_allq/" path5 <- "/Volumes/LaCie/demographic_simulation_paper/ms_and_vcf_simulation_output_20200504/" # specifying metric input files needed scenario_vec <- c("scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_1_hard_sweep.*.sample_p2.80K.fvec", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_1_hard_sweep.*.sample_p2.80K.fvec", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_1_hard_sweep.*.sample_p2.80K.fvec", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_1_hard_sweep.*.sample_p2.80K.fvec", "scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_1_soft_sweep.*.sample_p2.80K.fvec", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_1_soft_sweep.*.sample_p2.80K.fvec", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_1_soft_sweep.*.sample_p2.80K.fvec", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_1_soft_sweep.*.sample_p2.80K.fvec", "scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_01_hard_sweep.*.sample_p2.80K.fvec", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_01_hard_sweep.*.sample_p2.80K.fvec", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_01_hard_sweep.*.sample_p2.80K.fvec", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_01_hard_sweep.*.sample_p2.80K.fvec", "scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_01_soft_sweep.*.sample_p2.80K.fvec", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_01_soft_sweep.*.sample_p2.80K.fvec", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_01_soft_sweep.*.sample_p2.80K.fvec", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_01_soft_sweep.*.sample_p2.80K.fvec", "scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_001_hard_sweep.*.sample_p2.80K.fvec", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_001_hard_sweep.*.sample_p2.80K.fvec", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_001_hard_sweep.*.sample_p2.80K.fvec", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_001_hard_sweep.*.sample_p2.80K.fvec", "scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_001_soft_sweep.*.sample_p2.80K.fvec", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_001_soft_sweep.*.sample_p2.80K.fvec", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_001_soft_sweep.*.sample_p2.80K.fvec", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_001_soft_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_1_hard_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_1_hard_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_1_hard_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_1_hard_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_1_soft_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_1_soft_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_1_soft_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_1_soft_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_01_hard_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_01_hard_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_01_hard_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_01_hard_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_01_soft_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_01_soft_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_01_soft_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_01_soft_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_001_hard_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_001_hard_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_001_hard_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_001_hard_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_001_soft_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_001_soft_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_001_soft_sweep.*.sample_p2.80K.fvec", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_001_soft_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_1_hard_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_1_hard_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_1_hard_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_1_hard_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_1_soft_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_1_soft_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_1_soft_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_1_soft_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_01_hard_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_01_hard_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_01_hard_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_01_hard_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_01_soft_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_01_soft_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_01_soft_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_01_soft_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_001_hard_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_001_hard_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_001_hard_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_001_hard_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_001_soft_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_001_soft_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_001_soft_sweep.*.sample_p2.80K.fvec", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_001_soft_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_1_hard_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_1_hard_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_1_hard_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_1_hard_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_1_soft_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_1_soft_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_1_soft_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_1_soft_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_01_hard_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_01_hard_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_01_hard_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_01_hard_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_01_soft_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_01_soft_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_01_soft_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_01_soft_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_001_hard_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_001_hard_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_001_hard_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_001_hard_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_001_soft_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_001_soft_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_001_soft_sweep.*.sample_p2.80K.fvec", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_001_soft_sweep.*.sample_p2.80K.fvec") scenario_vec_p3 <- c("scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_1_hard_sweep.*.sample_p3.80K.fvec", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_1_hard_sweep.*.sample_p3.80K.fvec", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_1_hard_sweep.*.sample_p3.80K.fvec", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_1_hard_sweep.*.sample_p3.80K.fvec", "scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_1_soft_sweep.*.sample_p3.80K.fvec", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_1_soft_sweep.*.sample_p3.80K.fvec", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_1_soft_sweep.*.sample_p3.80K.fvec", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_1_soft_sweep.*.sample_p3.80K.fvec", "scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_01_hard_sweep.*.sample_p3.80K.fvec", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_01_hard_sweep.*.sample_p3.80K.fvec", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_01_hard_sweep.*.sample_p3.80K.fvec", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_01_hard_sweep.*.sample_p3.80K.fvec", "scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_01_soft_sweep.*.sample_p3.80K.fvec", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_01_soft_sweep.*.sample_p3.80K.fvec", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_01_soft_sweep.*.sample_p3.80K.fvec", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_01_soft_sweep.*.sample_p3.80K.fvec", "scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_001_hard_sweep.*.sample_p3.80K.fvec", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_001_hard_sweep.*.sample_p3.80K.fvec", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_001_hard_sweep.*.sample_p3.80K.fvec", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_001_hard_sweep.*.sample_p3.80K.fvec", "scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_001_soft_sweep.*.sample_p3.80K.fvec", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_001_soft_sweep.*.sample_p3.80K.fvec", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_001_soft_sweep.*.sample_p3.80K.fvec", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_001_soft_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_1_hard_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_1_hard_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_1_hard_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_1_hard_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_1_soft_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_1_soft_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_1_soft_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_1_soft_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_01_hard_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_01_hard_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_01_hard_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_01_hard_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_01_soft_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_01_soft_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_01_soft_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_01_soft_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_001_hard_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_001_hard_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_001_hard_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_001_hard_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_001_soft_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_001_soft_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_001_soft_sweep.*.sample_p3.80K.fvec", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_001_soft_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_1_hard_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_1_hard_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_1_hard_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_1_hard_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_1_soft_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_1_soft_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_1_soft_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_1_soft_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_01_hard_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_01_hard_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_01_hard_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_01_hard_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_01_soft_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_01_soft_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_01_soft_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_01_soft_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_001_hard_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_001_hard_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_001_hard_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_001_hard_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_001_soft_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_001_soft_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_001_soft_sweep.*.sample_p3.80K.fvec", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_001_soft_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_1_hard_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_1_hard_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_1_hard_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_1_hard_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_1_soft_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_1_soft_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_1_soft_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_1_soft_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_01_hard_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_01_hard_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_01_hard_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_01_hard_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_01_soft_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_01_soft_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_01_soft_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_01_soft_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_001_hard_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_001_hard_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_001_hard_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_001_hard_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_001_soft_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_001_soft_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_001_soft_sweep.*.sample_p3.80K.fvec", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_001_soft_sweep.*.sample_p3.80K.fvec") beta_vec <- c("scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_1_hard_sweep.*.sample.rds", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_1_hard_sweep.*.sample.rds", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_1_hard_sweep.*.sample.rds", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_1_hard_sweep.*.sample.rds", "scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_1_soft_sweep.*.sample.rds", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_1_soft_sweep.*.sample.rds", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_1_soft_sweep.*.sample.rds", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_1_soft_sweep.*.sample.rds", "scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_01_hard_sweep.*.sample.rds", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_01_hard_sweep.*.sample.rds", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_01_hard_sweep.*.sample.rds", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_01_hard_sweep.*.sample.rds", "scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_01_soft_sweep.*.sample.rds", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_01_soft_sweep.*.sample.rds", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_01_soft_sweep.*.sample.rds", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_01_soft_sweep.*.sample.rds", "scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_001_hard_sweep.*.sample.rds", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_001_hard_sweep.*.sample.rds", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_001_hard_sweep.*.sample.rds", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_001_hard_sweep.*.sample.rds", "scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_001_soft_sweep.*.sample.rds", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_001_soft_sweep.*.sample.rds", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_001_soft_sweep.*.sample.rds", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_001_soft_sweep.*.sample.rds", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_1_hard_sweep.*.sample.rds", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_1_hard_sweep.*.sample.rds", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_1_hard_sweep.*.sample.rds", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_1_hard_sweep.*.sample.rds", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_1_soft_sweep.*.sample.rds", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_1_soft_sweep.*.sample.rds", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_1_soft_sweep.*.sample.rds", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_1_soft_sweep.*.sample.rds", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_01_hard_sweep.*.sample.rds", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_01_hard_sweep.*.sample.rds", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_01_hard_sweep.*.sample.rds", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_01_hard_sweep.*.sample.rds", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_01_soft_sweep.*.sample.rds", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_01_soft_sweep.*.sample.rds", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_01_soft_sweep.*.sample.rds", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_01_soft_sweep.*.sample.rds", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_001_hard_sweep.*.sample.rds", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_001_hard_sweep.*.sample.rds", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_001_hard_sweep.*.sample.rds", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_001_hard_sweep.*.sample.rds", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_001_soft_sweep.*.sample.rds", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_001_soft_sweep.*.sample.rds", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_001_soft_sweep.*.sample.rds", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_001_soft_sweep.*.sample.rds", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_1_hard_sweep.*.sample.rds", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_1_hard_sweep.*.sample.rds", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_1_hard_sweep.*.sample.rds", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_1_hard_sweep.*.sample.rds", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_1_soft_sweep.*.sample.rds", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_1_soft_sweep.*.sample.rds", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_1_soft_sweep.*.sample.rds", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_1_soft_sweep.*.sample.rds", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_01_hard_sweep.*.sample.rds", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_01_hard_sweep.*.sample.rds", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_01_hard_sweep.*.sample.rds", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_01_hard_sweep.*.sample.rds", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_01_soft_sweep.*.sample.rds", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_01_soft_sweep.*.sample.rds", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_01_soft_sweep.*.sample.rds", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_01_soft_sweep.*.sample.rds", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_001_hard_sweep.*.sample.rds", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_001_hard_sweep.*.sample.rds", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_001_hard_sweep.*.sample.rds", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_001_hard_sweep.*.sample.rds", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_001_soft_sweep.*.sample.rds", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_001_soft_sweep.*.sample.rds", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_001_soft_sweep.*.sample.rds", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_001_soft_sweep.*.sample.rds", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_1_hard_sweep.*.sample.rds", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_1_hard_sweep.*.sample.rds", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_1_hard_sweep.*.sample.rds", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_1_hard_sweep.*.sample.rds", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_1_soft_sweep.*.sample.rds", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_1_soft_sweep.*.sample.rds", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_1_soft_sweep.*.sample.rds", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_1_soft_sweep.*.sample.rds", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_01_hard_sweep.*.sample.rds", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_01_hard_sweep.*.sample.rds", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_01_hard_sweep.*.sample.rds", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_01_hard_sweep.*.sample.rds", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_01_soft_sweep.*.sample.rds", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_01_soft_sweep.*.sample.rds", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_01_soft_sweep.*.sample.rds", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_01_soft_sweep.*.sample.rds", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_001_hard_sweep.*.sample.rds", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_001_hard_sweep.*.sample.rds", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_001_hard_sweep.*.sample.rds", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_001_hard_sweep.*.sample.rds", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_001_soft_sweep.*.sample.rds", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_001_soft_sweep.*.sample.rds", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_001_soft_sweep.*.sample.rds", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_001_soft_sweep.*.sample.rds") fst_vec <- c("scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_1_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_1_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_1_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_1_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_1_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_1_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_1_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_1_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_01_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_01_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_01_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_01_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_01_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_01_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_01_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_01_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_001_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_001_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_001_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_001_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_001_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario1_sim_sel0_1_constant_equal_popsize_migrate0_001_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_001_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_001_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_1_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_1_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_1_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_1_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_1_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_1_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_1_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_1_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_01_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_01_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_01_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_01_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_01_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_01_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_01_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_01_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_001_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_001_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_001_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_001_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_001_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_constant_equal_popsize_migrate0_001_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_001_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_001_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_1_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_1_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_1_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_1_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_1_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_1_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_1_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_1_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_01_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_01_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_01_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_01_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_01_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_01_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_01_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_01_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_001_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_001_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_001_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_001_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_001_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_constant_equal_popsize_migrate0_001_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_001_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_001_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_1_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_1_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_1_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_1_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_1_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_1_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_1_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_1_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_01_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_01_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_01_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_01_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_01_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_01_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_01_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_01_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_001_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_001_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_001_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_001_hard_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_001_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_constant_equal_popsize_migrate0_001_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_001_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst", "scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_001_soft_sweep.*.sample.80001.noNA.fst.windowed.weir.fst") scenario_vec_neutral <- c("scenario1_full_size_neutral*constant*0_5K_2K_*migrate0_1*.ms.sample_p2.80K.fvec", "scenario1_full_size_neutral*constant*0_5K_sel2_*migrate0_1*.ms.sample_p2.80K.fvec", "scenario1_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_1*.ms.sample_p2.80K.fvec", "scenario1_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_1*.ms.sample_p2.80K.fvec", "scenario1_full_size_neutral*constant*0_5K_2K_*migrate0_01*.ms.sample_p2.80K.fvec", "scenario1_full_size_neutral*constant*0_5K_sel2_*migrate0_01*.ms.sample_p2.80K.fvec", "scenario1_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_01*.ms.sample_p2.80K.fvec", "scenario1_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_01*.ms.sample_p2.80K.fvec", "scenario1_full_size_neutral*constant*0_5K_2K_*migrate0_001*.ms.sample_p2.80K.fvec", "scenario1_full_size_neutral*constant*0_5K_sel2_*migrate0_001*.ms.sample_p2.80K.fvec", "scenario1_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_001*.ms.sample_p2.80K.fvec", "scenario1_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_001*.ms.sample_p2.80K.fvec", "scenario2_full_size_neutral*constant*0_5K_2K_*migrate0_1*.ms.sample_p2.80K.fvec", "scenario2_full_size_neutral*constant*0_5K_sel2_*migrate0_1*.ms.sample_p2.80K.fvec", "scenario2_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_1*.ms.sample_p2.80K.fvec", "scenario2_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_1*.ms.sample_p2.80K.fvec", "scenario2_full_size_neutral*constant*0_5K_2K_*migrate0_01*.ms.sample_p2.80K.fvec", "scenario2_full_size_neutral*constant*0_5K_sel2_*migrate0_01*.ms.sample_p2.80K.fvec", "scenario2_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_01*.ms.sample_p2.80K.fvec", "scenario2_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_01*.ms.sample_p2.80K.fvec", "scenario2_full_size_neutral*constant*0_5K_2K_*migrate0_001*.ms.sample_p2.80K.fvec", "scenario2_full_size_neutral*constant*0_5K_sel2_*migrate0_001*.ms.sample_p2.80K.fvec", "scenario2_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_001*.ms.sample_p2.80K.fvec", "scenario2_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_001*.ms.sample_p2.80K.fvec", "scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_1*.ms.sample_p2.80K.fvec", "scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_1*.ms.sample_p2.80K.fvec", "scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_1*.ms.sample_p2.80K.fvec", "scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_1*.ms.sample_p2.80K.fvec", "scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_01*.ms.sample_p2.80K.fvec", "scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_01*.ms.sample_p2.80K.fvec", "scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_01*.ms.sample_p2.80K.fvec", "scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_01*.ms.sample_p2.80K.fvec", "scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_001*.ms.sample_p2.80K.fvec", "scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_001*.ms.sample_p2.80K.fvec", "scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_001*.ms.sample_p2.80K.fvec", "scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_001*.ms.sample_p2.80K.fvec", "diff_scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_1*.ms.sample_p2.80K.fvec", "diff_scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_1*.ms.sample_p2.80K.fvec", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_1*.ms.sample_p2.80K.fvec", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_1*.ms.sample_p2.80K.fvec", "diff_scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_01*.ms.sample_p2.80K.fvec", "diff_scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_01*.ms.sample_p2.80K.fvec", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_01*.ms.sample_p2.80K.fvec", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_01*.ms.sample_p2.80K.fvec", "diff_scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_001*.ms.sample_p2.80K.fvec", "diff_scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_001*.ms.sample_p2.80K.fvec", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_001*.ms.sample_p2.80K.fvec", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_001*.ms.sample_p2.80K.fvec") scenario_vec_p3_neutral <- c("scenario1_full_size_neutral*constant*0_5K_2K_*migrate0_1*.ms.sample_p3.80K.fvec", "scenario1_full_size_neutral*constant*0_5K_sel2_*migrate0_1*.ms.sample_p3.80K.fvec", "scenario1_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_1*.ms.sample_p3.80K.fvec", "scenario1_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_1*.ms.sample_p3.80K.fvec", "scenario1_full_size_neutral*constant*0_5K_2K_*migrate0_01*.ms.sample_p3.80K.fvec", "scenario1_full_size_neutral*constant*0_5K_sel2_*migrate0_01*.ms.sample_p3.80K.fvec", "scenario1_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_01*.ms.sample_p3.80K.fvec", "scenario1_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_01*.ms.sample_p3.80K.fvec", "scenario1_full_size_neutral*constant*0_5K_2K_*migrate0_001*.ms.sample_p3.80K.fvec", "scenario1_full_size_neutral*constant*0_5K_sel2_*migrate0_001*.ms.sample_p3.80K.fvec", "scenario1_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_001*.ms.sample_p3.80K.fvec", "scenario1_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_001*.ms.sample_p3.80K.fvec", "scenario2_full_size_neutral*constant*0_5K_2K_*migrate0_1*.ms.sample_p3.80K.fvec", "scenario2_full_size_neutral*constant*0_5K_sel2_*migrate0_1*.ms.sample_p3.80K.fvec", "scenario2_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_1*.ms.sample_p3.80K.fvec", "scenario2_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_1*.ms.sample_p3.80K.fvec", "scenario2_full_size_neutral*constant*0_5K_2K_*migrate0_01*.ms.sample_p3.80K.fvec", "scenario2_full_size_neutral*constant*0_5K_sel2_*migrate0_01*.ms.sample_p3.80K.fvec", "scenario2_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_01*.ms.sample_p3.80K.fvec", "scenario2_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_01*.ms.sample_p3.80K.fvec", "scenario2_full_size_neutral*constant*0_5K_2K_*migrate0_001*.ms.sample_p3.80K.fvec", "scenario2_full_size_neutral*constant*0_5K_sel2_*migrate0_001*.ms.sample_p3.80K.fvec", "scenario2_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_001*.ms.sample_p3.80K.fvec", "scenario2_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_001*.ms.sample_p3.80K.fvec", "scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_1*.ms.sample_p3.80K.fvec", "scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_1*.ms.sample_p3.80K.fvec", "scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_1*.ms.sample_p3.80K.fvec", "scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_1*.ms.sample_p3.80K.fvec", "scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_01*.ms.sample_p3.80K.fvec", "scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_01*.ms.sample_p3.80K.fvec", "scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_01*.ms.sample_p3.80K.fvec", "scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_01*.ms.sample_p3.80K.fvec", "scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_001*.ms.sample_p3.80K.fvec", "scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_001*.ms.sample_p3.80K.fvec", "scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_001*.ms.sample_p3.80K.fvec", "scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_001*.ms.sample_p3.80K.fvec", "diff_scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_1*.ms.sample_p3.80K.fvec", "diff_scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_1*.ms.sample_p3.80K.fvec", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_1*.ms.sample_p3.80K.fvec", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_1*.ms.sample_p3.80K.fvec", "diff_scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_01*.ms.sample_p3.80K.fvec", "diff_scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_01*.ms.sample_p3.80K.fvec", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_01*.ms.sample_p3.80K.fvec", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_01*.ms.sample_p3.80K.fvec", "diff_scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_001*.ms.sample_p3.80K.fvec", "diff_scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_001*.ms.sample_p3.80K.fvec", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_001*.ms.sample_p3.80K.fvec", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_001*.ms.sample_p3.80K.fvec") fst_vec_neutral <- c("scenario1_full_size_neutral*constant*0_5K_2K_*migrate0_1*.80001.noNA.fst.windowed.weir.fst", "scenario1_full_size_neutral*constant*0_5K_sel2_*migrate0_1*.80001.noNA.fst.windowed.weir.fst", "scenario1_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_1*.80001.noNA.fst.windowed.weir.fst", "scenario1_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_1*.80001.noNA.fst.windowed.weir.fst", "scenario1_full_size_neutral*constant*0_5K_2K_*migrate0_01*.80001.noNA.fst.windowed.weir.fst", "scenario1_full_size_neutral*constant*0_5K_sel2_*migrate0_01*.80001.noNA.fst.windowed.weir.fst", "scenario1_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_01*.80001.noNA.fst.windowed.weir.fst", "scenario1_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_01*.80001.noNA.fst.windowed.weir.fst", "scenario1_full_size_neutral*constant*0_5K_2K_*migrate0_001*.80001.noNA.fst.windowed.weir.fst", "scenario1_full_size_neutral*constant*0_5K_sel2_*migrate0_001*.80001.noNA.fst.windowed.weir.fst", "scenario1_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_001*.80001.noNA.fst.windowed.weir.fst", "scenario1_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_001*.80001.noNA.fst.windowed.weir.fst", "scenario2_full_size_neutral*constant*0_5K_2K_*migrate0_1*.80001.noNA.fst.windowed.weir.fst", "scenario2_full_size_neutral*constant*0_5K_sel2_*migrate0_1*.80001.noNA.fst.windowed.weir.fst", "scenario2_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_1*.80001.noNA.fst.windowed.weir.fst", "scenario2_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_1*.80001.noNA.fst.windowed.weir.fst", "scenario2_full_size_neutral*constant*0_5K_2K_*migrate0_01*.80001.noNA.fst.windowed.weir.fst", "scenario2_full_size_neutral*constant*0_5K_sel2_*migrate0_01*.80001.noNA.fst.windowed.weir.fst", "scenario2_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_01*.80001.noNA.fst.windowed.weir.fst", "scenario2_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_01*.80001.noNA.fst.windowed.weir.fst", "scenario2_full_size_neutral*constant*0_5K_2K_*migrate0_001*.80001.noNA.fst.windowed.weir.fst", "scenario2_full_size_neutral*constant*0_5K_sel2_*migrate0_001*.80001.noNA.fst.windowed.weir.fst", "scenario2_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_001*.80001.noNA.fst.windowed.weir.fst", "scenario2_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_001*.80001.noNA.fst.windowed.weir.fst", "scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_1*.80001.noNA.fst.windowed.weir.fst", "scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_1*.80001.noNA.fst.windowed.weir.fst", "scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_1*.80001.noNA.fst.windowed.weir.fst", "scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_1*.80001.noNA.fst.windowed.weir.fst", "scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_01*.80001.noNA.fst.windowed.weir.fst", "scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_01*.80001.noNA.fst.windowed.weir.fst", "scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_01*.80001.noNA.fst.windowed.weir.fst", "scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_01*.80001.noNA.fst.windowed.weir.fst", "scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_001*.80001.noNA.fst.windowed.weir.fst", "scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_001*.80001.noNA.fst.windowed.weir.fst", "scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_001*.80001.noNA.fst.windowed.weir.fst", "scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_001*.80001.noNA.fst.windowed.weir.fst", "diff_scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_1*.80001.noNA.fst.windowed.weir.fst", "diff_scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_1*.80001.noNA.fst.windowed.weir.fst", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_1*.80001.noNA.fst.windowed.weir.fst", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_1*.80001.noNA.fst.windowed.weir.fst", "diff_scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_01*.80001.noNA.fst.windowed.weir.fst", "diff_scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_01*.80001.noNA.fst.windowed.weir.fst", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_01*.80001.noNA.fst.windowed.weir.fst", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_01*.80001.noNA.fst.windowed.weir.fst", "diff_scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_001*.80001.noNA.fst.windowed.weir.fst", "diff_scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_001*.80001.noNA.fst.windowed.weir.fst", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_001*.80001.noNA.fst.windowed.weir.fst", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_001*.80001.noNA.fst.windowed.weir.fst") beta_vec_neutral <- c("scenario1_full_size_neutral*constant*0_5K_2K_*migrate0_1*.ms.sample.comb.rds", "scenario1_full_size_neutral*constant*0_5K_sel2_*migrate0_1*.ms.sample.comb.rds", "scenario1_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_1*.ms.sample.comb.rds", "scenario1_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_1*.ms.sample.comb.rds", "scenario1_full_size_neutral*constant*0_5K_2K_*migrate0_01*.ms.sample.comb.rds", "scenario1_full_size_neutral*constant*0_5K_sel2_*migrate0_01*.ms.sample.comb.rds", "scenario1_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_01*.ms.sample.comb.rds", "scenario1_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_01*.ms.sample.comb.rds", "scenario1_full_size_neutral*constant*0_5K_2K_*migrate0_001*.ms.sample.comb.rds", "scenario1_full_size_neutral*constant*0_5K_sel2_*migrate0_001*.ms.sample.comb.rds", "scenario1_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_001*.ms.sample.comb.rds", "scenario1_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_001*.ms.sample.comb.rds", "scenario2_full_size_neutral*constant*0_5K_2K_*migrate0_1*.ms.sample.comb.rds", "scenario2_full_size_neutral*constant*0_5K_sel2_*migrate0_1*.ms.sample.comb.rds", "scenario2_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_1*.ms.sample.comb.rds", "scenario2_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_1*.ms.sample.comb.rds", "scenario2_full_size_neutral*constant*0_5K_2K_*migrate0_01*.ms.sample.comb.rds", "scenario2_full_size_neutral*constant*0_5K_sel2_*migrate0_01*.ms.sample.comb.rds", "scenario2_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_01*.ms.sample.comb.rds", "scenario2_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_01*.ms.sample.comb.rds", "scenario2_full_size_neutral*constant*0_5K_2K_*migrate0_001*.ms.sample.comb.rds", "scenario2_full_size_neutral*constant*0_5K_sel2_*migrate0_001*.ms.sample.comb.rds", "scenario2_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_001*.ms.sample.comb.rds", "scenario2_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_001*.ms.sample.comb.rds", "scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_1*.ms.sample.comb.rds", "scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_1*.ms.sample.comb.rds", "scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_1*.ms.sample.comb.rds", "scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_1*.ms.sample.comb.rds", "scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_01*.ms.sample.comb.rds", "scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_01*.ms.sample.comb.rds", "scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_01*.ms.sample.comb.rds", "scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_01*.ms.sample.comb.rds", "scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_001*.ms.sample.comb.rds", "scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_001*.ms.sample.comb.rds", "scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_001*.ms.sample.comb.rds", "scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_001*.ms.sample.comb.rds", "diff_scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_1*.ms.sample.comb.rds", "diff_scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_1*.ms.sample.comb.rds", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_1*.ms.sample.comb.rds", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_1*.ms.sample.comb.rds", "diff_scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_01*.ms.sample.comb.rds", "diff_scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_01*.ms.sample.comb.rds", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_01*.ms.sample.comb.rds", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_01*.ms.sample.comb.rds", "diff_scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_001*.ms.sample.comb.rds", "diff_scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_001*.ms.sample.comb.rds", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_001*.ms.sample.comb.rds", "diff_scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_001*.ms.sample.comb.rds") hscan_vec <- c("hscan_scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_1_hard_sweep.*.sample_p2.ms", "hscan_scenario1_sim_sel0_1_constant_equal_popsize_migrate0_1_hard_sweep.*.sample_p2.ms", "hscan_scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_1_hard_sweep.*.sample_p2.ms", "hscan_scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_1_hard_sweep.*.sample_p2.ms", "hscan_scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_1_soft_sweep.*.sample_p2.ms", "hscan_scenario1_sim_sel0_1_constant_equal_popsize_migrate0_1_soft_sweep.*.sample_p2.ms", "hscan_scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_1_soft_sweep.*.sample_p2.ms", "hscan_scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_1_soft_sweep.*.sample_p2.ms", "hscan_scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_01_hard_sweep.*.sample_p2.ms", "hscan_scenario1_sim_sel0_1_constant_equal_popsize_migrate0_01_hard_sweep.*.sample_p2.ms", "hscan_scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_01_hard_sweep.*.sample_p2.ms", "hscan_scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_01_hard_sweep.*.sample_p2.ms", "hscan_scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_01_soft_sweep.*.sample_p2.ms", "hscan_scenario1_sim_sel0_1_constant_equal_popsize_migrate0_01_soft_sweep.*.sample_p2.ms", "hscan_scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_01_soft_sweep.*.sample_p2.ms", "hscan_scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_01_soft_sweep.*.sample_p2.ms", "hscan_scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_001_hard_sweep.*.sample_p2.ms", "hscan_scenario1_sim_sel0_1_constant_equal_popsize_migrate0_001_hard_sweep.*.sample_p2.ms", "hscan_scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_001_hard_sweep.*.sample_p2.ms", "hscan_scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_001_hard_sweep.*.sample_p2.ms", "hscan_scenario1_sim_sel0_1_constant_unequal_popsize_migrate0_001_soft_sweep.*.sample_p2.ms", "hscan_scenario1_sim_sel0_1_constant_equal_popsize_migrate0_001_soft_sweep.*.sample_p2.ms", "hscan_scenario1_sim_sel0_1_increasing_unequal_popsize_migrate0_001_soft_sweep.*.sample_p2.ms", "hscan_scenario1_sim_sel0_1_increasing_equal_popsize_migrate0_001_soft_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_1_hard_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_constant_equal_popsize_migrate0_1_hard_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_1_hard_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_1_hard_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_1_soft_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_constant_equal_popsize_migrate0_1_soft_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_1_soft_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_1_soft_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_01_hard_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_constant_equal_popsize_migrate0_01_hard_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_01_hard_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_01_hard_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_01_soft_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_constant_equal_popsize_migrate0_01_soft_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_01_soft_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_01_soft_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_001_hard_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_constant_equal_popsize_migrate0_001_hard_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_001_hard_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_001_hard_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_constant_unequal_popsize_migrate0_001_soft_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_constant_equal_popsize_migrate0_001_soft_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_increasing_unequal_popsize_migrate0_001_soft_sweep.*.sample_p2.ms", "hscan_scenario2_sim_sel0_1_increasing_equal_popsize_migrate0_001_soft_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_1_hard_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_constant_equal_popsize_migrate0_1_hard_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_1_hard_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_1_hard_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_1_soft_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_constant_equal_popsize_migrate0_1_soft_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_1_soft_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_1_soft_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_01_hard_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_constant_equal_popsize_migrate0_01_hard_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_01_hard_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_01_hard_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_01_soft_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_constant_equal_popsize_migrate0_01_soft_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_01_soft_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_01_soft_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_001_hard_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_constant_equal_popsize_migrate0_001_hard_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_001_hard_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_001_hard_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_constant_unequal_popsize_migrate0_001_soft_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_constant_equal_popsize_migrate0_001_soft_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_increasing_unequal_popsize_migrate0_001_soft_sweep.*.sample_p2.ms", "hscan_scenario3_sim_sel0_1_increasing_equal_popsize_migrate0_001_soft_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_1_hard_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_constant_equal_popsize_migrate0_1_hard_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_1_hard_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_1_hard_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_1_soft_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_constant_equal_popsize_migrate0_1_soft_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_1_soft_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_1_soft_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_01_hard_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_constant_equal_popsize_migrate0_01_hard_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_01_hard_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_01_hard_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_01_soft_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_constant_equal_popsize_migrate0_01_soft_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_01_soft_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_01_soft_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_001_hard_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_constant_equal_popsize_migrate0_001_hard_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_001_hard_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_001_hard_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_constant_unequal_popsize_migrate0_001_soft_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_constant_equal_popsize_migrate0_001_soft_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_increasing_unequal_popsize_migrate0_001_soft_sweep.*.sample_p2.ms", "hscan_scenario4_sim_sel0_1_increasing_equal_popsize_migrate0_001_soft_sweep.*.sample_p2.ms") hscan_vec_neutral <- c("hscan_scenario1_full_size_neutral*constant*0_5K_2K_*migrate0_1*.ms.sample_p2.txt", "hscan_scenario1_full_size_neutral*constant*0_5K_sel2_*migrate0_1*.ms.sample_p2.txt", "hscan_scenario1_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_1*.ms.sample_p2.txt", "hscan_scenario1_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_1*.ms.sample_p2.txt", "hscan_scenario1_full_size_neutral*constant*0_5K_2K_*migrate0_01*.ms.sample_p2.txt", "hscan_scenario1_full_size_neutral*constant*0_5K_sel2_*migrate0_01*.ms.sample_p2.txt", "hscan_scenario1_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_01*.ms.sample_p2.txt", "hscan_scenario1_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_01*.ms.sample_p2.txt", "hscan_scenario1_full_size_neutral*constant*0_5K_2K_*migrate0_001*.ms.sample_p2.txt", "hscan_scenario1_full_size_neutral*constant*0_5K_sel2_*migrate0_001*.ms.sample_p2.txt", "hscan_scenario1_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_001*.ms.sample_p2.txt", "hscan_scenario1_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_001*.ms.sample_p2.txt", "hscan_scenario2_full_size_neutral*constant*0_5K_2K_*migrate0_1*.ms.sample_p2.txt", "hscan_scenario2_full_size_neutral*constant*0_5K_sel2_*migrate0_1*.ms.sample_p2.txt", "hscan_scenario2_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_1*.ms.sample_p2.txt", "hscan_scenario2_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_1*.ms.sample_p2.txt", "hscan_scenario2_full_size_neutral*constant*0_5K_2K_*migrate0_01*.ms.sample_p2.txt", "hscan_scenario2_full_size_neutral*constant*0_5K_sel2_*migrate0_01*.ms.sample_p2.txt", "hscan_scenario2_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_01*.ms.sample_p2.txt", "hscan_scenario2_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_01*.ms.sample_p2.txt", "hscan_scenario2_full_size_neutral*constant*0_5K_2K_*migrate0_001*.ms.sample_p2.txt", "hscan_scenario2_full_size_neutral*constant*0_5K_sel2_*migrate0_001*.ms.sample_p2.txt", "hscan_scenario2_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_001*.ms.sample_p2.txt", "hscan_scenario2_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_001*.ms.sample_p2.txt", "hscan_scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_1*.ms.sample_p2.txt", "hscan_scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_1*.ms.sample_p2.txt", "hscan_scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_1*.ms.sample_p2.txt", "hscan_scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_1*.ms.sample_p2.txt", "hscan_scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_01*.ms.sample_p2.txt", "hscan_scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_01*.ms.sample_p2.txt", "hscan_scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_01*.ms.sample_p2.txt", "hscan_scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_01*.ms.sample_p2.txt", "hscan_scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_001*.ms.sample_p2.txt", "hscan_scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_001*.ms.sample_p2.txt", "hscan_scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_001*.ms.sample_p2.txt", "hscan_scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_001*.ms.sample_p2.txt", "hscan_diff_scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_1*.ms.sample_p2.txt", "hscan_diff_scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_1*.ms.sample_p2.txt", "hscan_diff_scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_1*.ms.sample_p2.txt", "hscan_diff_scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_1*.ms.sample_p2.txt", "hscan_diff_scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_01*.ms.sample_p2.txt", "hscan_diff_scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_01*.ms.sample_p2.txt", "hscan_diff_scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_01*.ms.sample_p2.txt", "hscan_diff_scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_01*.ms.sample_p2.txt", "hscan_diff_scenario3_full_size_neutral*constant*0_5K_2K_*migrate0_001*.ms.sample_p2.txt", "hscan_diff_scenario3_full_size_neutral*constant*0_5K_sel2_*migrate0_001*.ms.sample_p2.txt", "hscan_diff_scenario3_full_size_neutral*bigchr_popsize0_5K_2K_*migrate0_001*.ms.sample_p2.txt", "hscan_diff_scenario3_full_size_neutral*bigchr_popsize0_5K_sel2_*migrate0_001*.ms.sample_p2.txt") title_vec <- c("scenario 1 hard sweep (constant, unequal popsize; migrate 0.1)", "scenario 1 hard sweep (constant, equal popsize; migrate 0.1)", "scenario 1 hard sweep (growing, unequal popsize; migrate 0.1)", "scenario 1 hard sweep (growing, equal popsize; migrate 0.1)", "scenario 1 soft sweep (constant, unequal popsize; migrate 0.1)", "scenario 1 soft sweep (constant, equal popsize; migrate 0.1)", "scenario 1 soft sweep (growing, unequal popsize; migrate 0.1)", "scenario 1 soft sweep (growing, equal popsize; migrate 0.1)", "scenario 1 hard sweep (constant, unequal popsize; migrate 0.01)", "scenario 1 hard sweep (constant, equal popsize; migrate 0.01)", "scenario 1 hard sweep (growing, unequal popsize; migrate 0.01)", "scenario 1 hard sweep (growing, equal popsize; migrate 0.01)", "scenario 1 soft sweep (constant, unequal popsize; migrate 0.01)", "scenario 1 soft sweep (constant, equal popsize; migrate 0.01)", "scenario 1 soft sweep (growing, unequal popsize; migrate 0.01)", "scenario 1 soft sweep (growing, equal popsize; migrate 0.01)", "scenario 1 hard sweep (constant, unequal popsize; migrate 0.001)", "scenario 1 hard sweep (constant, equal popsize; migrate 0.001)", "scenario 1 hard sweep (growing, unequal popsize; migrate 0.001)", "scenario 1 hard sweep (growing, equal popsize; migrate 0.001)", "scenario 1 soft sweep (constant, unequal popsize; migrate 0.001)", "scenario 1 soft sweep (constant, equal popsize; migrate 0.001)", "scenario 1 soft sweep (growing, unequal popsize; migrate 0.001)", "scenario 1 soft sweep (growing, equal popsize; migrate 0.001)", "scenario 2 hard sweep (constant, unequal popsize; migrate 0.1)", "scenario 2 hard sweep (constant, equal popsize; migrate 0.1)", "scenario 2 hard sweep (growing, unequal popsize; migrate 0.1)", "scenario 2 hard sweep (growing, equal popsize; migrate 0.1)", "scenario 2 soft sweep (constant, unequal popsize; migrate 0.1)", "scenario 2 soft sweep (constant, equal popsize; migrate 0.1)", "scenario 2 soft sweep (growing, unequal popsize; migrate 0.1)", "scenario 2 soft sweep (growing, equal popsize; migrate 0.1)", "scenario 2 hard sweep (constant, unequal popsize; migrate 0.01)", "scenario 2 hard sweep (constant, equal popsize; migrate 0.01)", "scenario 2 hard sweep (growing, unequal popsize; migrate 0.01)", "scenario 2 hard sweep (growing, equal popsize; migrate 0.01)", "scenario 2 soft sweep (constant, unequal popsize; migrate 0.01)", "scenario 2 soft sweep (constant, equal popsize; migrate 0.01)", "scenario 2 soft sweep (growing, unequal popsize; migrate 0.01)", "scenario 2 soft sweep (growing, equal popsize; migrate 0.01)", "scenario 2 hard sweep (constant, unequal popsize; migrate 0.001)", "scenario 2 hard sweep (constant, equal popsize; migrate 0.001)", "scenario 2 hard sweep (growing, unequal popsize; migrate 0.001)", "scenario 2 hard sweep (growing, equal popsize; migrate 0.001)", "scenario 2 soft sweep (constant, unequal popsize; migrate 0.001)", "scenario 2 soft sweep (constant, equal popsize; migrate 0.001)", "scenario 2 soft sweep (growing, unequal popsize; migrate 0.001)", "scenario 2 soft sweep (growing, equal popsize; migrate 0.001)", "scenario 3 hard sweep (constant, unequal popsize; migrate 0.1)", "scenario 3 hard sweep (constant, equal popsize; migrate 0.1)", "scenario 3 hard sweep (growing, unequal popsize; migrate 0.1)", "scenario 3 hard sweep (growing, equal popsize; migrate 0.1)", "scenario 3 soft sweep (constant, unequal popsize; migrate 0.1)", "scenario 3 soft sweep (constant, equal popsize; migrate 0.1)", "scenario 3 soft sweep (growing, unequal popsize; migrate 0.1)", "scenario 3 soft sweep (growing, equal popsize; migrate 0.1)", "scenario 3 hard sweep (constant, unequal popsize; migrate 0.01)", "scenario 3 hard sweep (constant, equal popsize; migrate 0.01)", "scenario 3 hard sweep (growing, unequal popsize; migrate 0.01)", "scenario 3 hard sweep (growing, equal popsize; migrate 0.01)", "scenario 3 soft sweep (constant, unequal popsize; migrate 0.01)", "scenario 3 soft sweep (constant, equal popsize; migrate 0.01)", "scenario 3 soft sweep (growing, unequal popsize; migrate 0.01)", "scenario 3 soft sweep (growing, equal popsize; migrate 0.01)", "scenario 3 hard sweep (constant, unequal popsize; migrate 0.001)", "scenario 3 hard sweep (constant, equal popsize; migrate 0.001)", "scenario 3 hard sweep (growing, unequal popsize; migrate 0.001)", "scenario 3 hard sweep (growing, equal popsize; migrate 0.001)", "scenario 3 soft sweep (constant, unequal popsize; migrate 0.001)", "scenario 3 soft sweep (constant, equal popsize; migrate 0.001)", "scenario 3 soft sweep (growing, unequal popsize; migrate 0.001)", "scenario 3 soft sweep (growing, equal popsize; migrate 0.001)", "scenario 4 hard sweep (constant, unequal popsize; migrate 0.1)", "scenario 4 hard sweep (constant, equal popsize; migrate 0.1)", "scenario 4 hard sweep (growing, unequal popsize; migrate 0.1)", "scenario 4 hard sweep (growing, equal popsize; migrate 0.1)", "scenario 4 soft sweep (constant, unequal popsize; migrate 0.1)", "scenario 4 soft sweep (constant, equal popsize; migrate 0.1)", "scenario 4 soft sweep (growing, unequal popsize; migrate 0.1)", "scenario 4 soft sweep (growing, equal popsize; migrate 0.1)", "scenario 4 hard sweep (constant, unequal popsize; migrate 0.01)", "scenario 4 hard sweep (constant, equal popsize; migrate 0.01)", "scenario 4 hard sweep (growing, unequal popsize; migrate 0.01)", "scenario 4 hard sweep (growing, equal popsize; migrate 0.01)", "scenario 4 soft sweep (constant, unequal popsize; migrate 0.01)", "scenario 4 soft sweep (constant, equal popsize; migrate 0.01)", "scenario 4 soft sweep (growing, unequal popsize; migrate 0.01)", "scenario 4 soft sweep (growing, equal popsize; migrate 0.01)", "scenario 4 hard sweep (constant, unequal popsize; migrate 0.001)", "scenario 4 hard sweep (constant, equal popsize; migrate 0.001)", "scenario 4 hard sweep (growing, unequal popsize; migrate 0.001)", "scenario 4 hard sweep (growing, equal popsize; migrate 0.001)", "scenario 4 soft sweep (constant, unequal popsize; migrate 0.001)", "scenario 4 soft sweep (constant, equal popsize; migrate 0.001)", "scenario 4 soft sweep (growing, unequal popsize; migrate 0.001)", "scenario 4 soft sweep (growing, equal popsize; migrate 0.001)") title_vec_neutral <- c("scenario 1 neutral (constant, unequal popsize; migrate 0.1)", "scenario 1 neutral (constant, equal popsize; migrate 0.1)", "scenario 1 neutral (growing, unequal popsize; migrate 0.1)", "scenario 1 neutral (growing, equal popsize; migrate 0.1)", "scenario 1 neutral (constant, unequal popsize; migrate 0.01)", "scenario 1 neutral (constant, equal popsize; migrate 0.01)", "scenario 1 neutral (growing, unequal popsize; migrate 0.01)", "scenario 1 neutral (growing, equal popsize; migrate 0.01)", "scenario 1 neutral (constant, unequal popsize; migrate 0.001)", "scenario 1 neutral (constant, equal popsize; migrate 0.001)", "scenario 1 neutral (growing, unequal popsize; migrate 0.001)", "scenario 1 neutral (growing, equal popsize; migrate 0.001)", "scenario 2 neutral (constant, unequal popsize; migrate 0.1)", "scenario 2 neutral (constant, equal popsize; migrate 0.1)", "scenario 2 neutral (growing, unequal popsize; migrate 0.1)", "scenario 2 neutral (growing, equal popsize; migrate 0.1)", "scenario 2 neutral (constant, unequal popsize; migrate 0.01)", "scenario 2 neutral (constant, equal popsize; migrate 0.01)", "scenario 2 neutral (growing, unequal popsize; migrate 0.01)", "scenario 2 neutral (growing, equal popsize; migrate 0.01)", "scenario 2 neutral (constant, unequal popsize; migrate 0.001)", "scenario 2 neutral (constant, equal popsize; migrate 0.001)", "scenario 2 neutral (growing, unequal popsize; migrate 0.001)", "scenario 2 neutral (growing, equal popsize; migrate 0.001)", "scenario 3 neutral (constant, unequal popsize; migrate 0.1)", "scenario 3 neutral (constant, equal popsize; migrate 0.1)", "scenario 3 neutral (growing, unequal popsize; migrate 0.1)", "scenario 3 neutral (growing, equal popsize; migrate 0.1)", "scenario 3 neutral (constant, unequal popsize; migrate 0.01)", "scenario 3 neutral (constant, equal popsize; migrate 0.01)", "scenario 3 neutral (growing, unequal popsize; migrate 0.01)", "scenario 3 neutral (growing, equal popsize; migrate 0.01)", "scenario 3 neutral (constant, unequal popsize; migrate 0.001)", "scenario 3 neutral (constant, equal popsize; migrate 0.001)", "scenario 3 neutral (growing, unequal popsize; migrate 0.001)", "scenario 3 neutral (growing, equal popsize; migrate 0.001)", "scenario 4 neutral (constant, unequal popsize; migrate 0.1)", "scenario 4 neutral (constant, equal popsize; migrate 0.1)", "scenario 4 neutral (growing, unequal popsize; migrate 0.1)", "scenario 4 neutral (growing, equal popsize; migrate 0.1)", "scenario 4 neutral (constant, unequal popsize; migrate 0.01)", "scenario 4 neutral (constant, equal popsize; migrate 0.01)", "scenario 4 neutral (growing, unequal popsize; migrate 0.01)", "scenario 4 neutral (growing, equal popsize; migrate 0.01)", "scenario 4 neutral (constant, unequal popsize; migrate 0.001)", "scenario 4 neutral (constant, equal popsize; migrate 0.001)", "scenario 4 neutral (growing, unequal popsize; migrate 0.001)", "scenario 4 neutral (growing, equal popsize; migrate 0.001)") ## performing a sliding window analysis based on the H-Scan input # function for sliding window analysis sliding_window <- function (input_frame) { input_frame$x <- (input_frame$x / 1000000) * 4000050 window_frame <- data.frame("pos" = NA, "hscan" = NA) for (j in seq(1, 4000050, 80001)) { hscan_mean <- NULL if ((j + 80000) > 4000050) { win_length <- 4000050 - j } else { win_length <- 80001 } # calculating the mean of the metric hscan_mean <- mean(input_frame$H[which(input_frame$x >= j & input_frame$x <= j + 80000)], na.rm = TRUE) window_frame$pos[nrow(window_frame)] <- j window_frame$hscan[nrow(window_frame)] <- hscan_mean window_frame <- rbind(window_frame, NA) } window_frame <- window_frame[-nrow(window_frame), ] return(window_frame) } for (i in 1:length(hscan_vec)) { glob_vec_pre1 <- Sys.glob(file.path(paste(path, hscan_vec[i], sep = ""))) for (j in 1:length(glob_vec_pre1)) { hscan_frame <- read.table(glob_vec_pre1[j], header = TRUE) hscan_window_frame <- sliding_window(hscan_frame) write.table(hscan_window_frame, paste(glob_vec_pre1[j], ".windows", sep = ""), quote = FALSE, sep = "\t", row.names = FALSE) } } for (i in 1:length(hscan_vec_neutral)) { glob_vec_pre2 <- Sys.glob(file.path(paste(path, hscan_vec_neutral[i], sep = ""))) for (j in 1:length(glob_vec_pre2)) { hscan_neutral_frame <- read.table(glob_vec_pre2[j], header = TRUE) hscan_neutral_window_frame <- sliding_window(hscan_neutral_frame) write.table(hscan_neutral_window_frame, paste(glob_vec_pre2[j], ".windows", sep = ""), quote = FALSE, sep = "\t", row.names = FALSE) } } ### # parsing and changing format of scenario names filename_vec <- sapply(title_vec, function(x) gsub(" ", "_", x)) filename_vec <- sapply(filename_vec, function(x) str_remove(x, "\\(")) filename_vec <- sapply(filename_vec, function(x) str_remove(x, "\\)")) filename_vec <- sapply(filename_vec, function(x) str_remove(x, ",")) filename_vec <- sapply(filename_vec, function(x) str_remove(x, ";")) filename_vec <- sapply(filename_vec, function(x) gsub("0\\.", "0_", x)) filename_vec_neutral <- sapply(title_vec_neutral, function(x) gsub(" ", "_", x)) filename_vec_neutral <- sapply(filename_vec_neutral, function(x) str_remove(x, "\\(")) filename_vec_neutral <- sapply(filename_vec_neutral, function(x) str_remove(x, "\\)")) filename_vec_neutral <- sapply(filename_vec_neutral, function(x) str_remove(x, ",")) filename_vec_neutral <- sapply(filename_vec_neutral, function(x) str_remove(x, ";")) filename_vec_neutral <- sapply(filename_vec_neutral, function(x) gsub("0\\.", "0_", x)) # initialising data frames maxp1_df <- as.data.frame(matrix(rep(0, 18 * length(title_vec)), length(title_vec), 18)) maxp2_df <- as.data.frame(matrix(rep(0, 18 * length(title_vec)), length(title_vec), 18)) minp1_df <- as.data.frame(matrix(rep(0, 18 * length(title_vec)), length(title_vec), 18)) minp2_df <- as.data.frame(matrix(rep(0, 18 * length(title_vec)), length(title_vec), 18)) for (i in 1:length(title_vec)) { print(i) glob_vec <- Sys.glob(file.path(paste(scenario_vec[i]))) glob_vec_p3 <- Sys.glob(file.path(paste(scenario_vec_p3[i]))) glob_vec2 <- Sys.glob(file.path(paste(path, fst_vec[i], sep = ""))) glob_vec3 <- Sys.glob(file.path(paste(path2, beta_vec[i], sep = ""))) glob_vec4 <- Sys.glob(file.path(paste(path, hscan_vec[i], ".windows", sep = ""))) df_list <- list() for (j in 1:100) { # reading in metric input files fv <- read.table(glob_vec[j], header = TRUE) fv_p3 <- read.table(glob_vec_p3[j], header = TRUE) statvars <- colsplit(colnames(fv), "_win", c("metric", "window")) fv <- unlist(sapply(fv, unlist)) fv_p3 <- unlist(sapply(fv_p3, unlist)) fv <- as.data.frame(matrix(fv, nrow = 50)) fv_p3 <- as.data.frame(matrix(fv_p3, nrow = 50)) rownames(fv) <- 1:50 rownames(fv_p3) <- 1:50 colnames(fv) <- unique(statvars[, 1]) colnames(fv_p3) <- unique(statvars[, 1]) hscan_frame <- read.table(glob_vec4[j], header = TRUE) fv <- data.frame(fv, "hscan" = hscan_frame[1:50, 2]) fv <- data.frame(fv, "delta_pi" = fv[1:50, 1] - fv_p3[1:50, 1]) fst_frame <- read.table(glob_vec2[j], header = TRUE) fv <- data.frame(fv, "Fst" = fst_frame[1:50, 6]) beta_frame <- readRDS(glob_vec3[j]) rownames(beta_frame) <- 1:50 fv <- data.frame(fv, "beta_q0" = beta_frame[1:50, 2]) fv <- data.frame(fv, "beta_q1" = beta_frame[1:50, 3]) fv <- data.frame(fv, "beta_q2" = beta_frame[1:50, 4]) fv <- data.frame(fv, "window_num" = (1:50) * 4000050 / 50) fv1 <- fv[1:25, ] fv2 <- fv[26:50, ] for (k in 1:(ncol(fv) - 1)) { part1 <- fv1[, k] part2 <- fv2[, k] maxp1_1 <- order(part1, decreasing = TRUE)[1] maxp1_2 <- order(part1, decreasing = TRUE)[2] maxp2_1 <- order(part2, decreasing = TRUE)[1] maxp2_2 <- order(part2, decreasing = TRUE)[2] minp1_1 <- order(part1, decreasing = FALSE)[1] minp1_2 <- order(part1, decreasing = FALSE)[2] minp2_1 <- order(part2, decreasing = FALSE)[1] minp2_2 <- order(part2, decreasing = FALSE)[2] # determining maxima and minima of metrics if (maxp1_1 == 13) { maxp1_df[i, k] <- maxp1_df[i, k] + 1 } if (maxp2_1 == 13) { maxp2_df[i, k] <- maxp2_df[i, k] + 1 } if (minp1_1 == 13) { minp1_df[i, k] <- minp1_df[i, k] + 1 } if (minp2_1 == 13) { minp2_df[i, k] <- minp2_df[i, k] + 1 } } } maxp1_df[i, ] <- maxp1_df[i, ] / length(glob_vec) maxp2_df[i, ] <- maxp2_df[i, ] / length(glob_vec) minp1_df[i, ] <- minp1_df[i, ] / length(glob_vec) minp2_df[i, ] <- minp2_df[i, ] / length(glob_vec) } colnames(maxp1_df) <- colnames(fv)[1:18] colnames(maxp2_df) <- colnames(fv)[1:18] colnames(minp1_df) <- colnames(fv)[1:18] colnames(minp2_df) <- colnames(fv)[1:18] maxp1_df_new <- data.frame("scenario" = title_vec, maxp1_df) maxp2_df_new <- data.frame("scenario" = title_vec, maxp2_df) minp1_df_new <- data.frame("scenario" = title_vec, minp1_df) minp2_df_new <- data.frame("scenario" = title_vec, minp2_df) write.table(maxp1_df_new, "rev_18_prediction_freq_maxp1_1in1_v2.tsv", quote = FALSE, sep = "\t", row.names = FALSE) write.table(maxp2_df_new, "rev_18_prediction_freq_maxp2_1in1_v2.tsv", quote = FALSE, sep = "\t", row.names = FALSE) write.table(minp1_df_new, "rev_18_prediction_freq_minp1_1in1_v2.tsv", quote = FALSE, sep = "\t", row.names = FALSE) write.table(minp2_df_new, "rev_18_prediction_freq_minp2_1in1_v2.tsv", quote = FALSE, sep = "\t", row.names = FALSE) maxp1_df_new <- data.frame(maxp1_df_new, "recombination" = "high") maxp2_df_new <- data.frame(maxp2_df_new, "recombination" = "low") minp1_df_new <- data.frame(minp1_df_new, "recombination" = "high") minp2_df_new <- data.frame(minp2_df_new, "recombination" = "low") max_df <- rbind(maxp1_df_new, maxp2_df_new) min_df <- rbind(minp1_df_new, minp2_df_new) rownames(max_df) <- 1:nrow(max_df) rownames(min_df) <- 1:nrow(min_df) max_df <- data.frame(max_df, "hard_sweep" = str_detect(max_df$scenario, "hard", negate = FALSE)) min_df <- data.frame(min_df, "hard_sweep" = str_detect(min_df$scenario, "hard", negate = FALSE)) # calculating means and standard deviations of maxima and minima proportions mean_max_df <- aggregate(max_df[, 2:19], list(max_df$recombination, max_df$hard_sweep), mean) mean_min_df <- aggregate(min_df[, 2:19], list(min_df$recombination, min_df$hard_sweep), mean) sd_max_df <- aggregate(max_df[, 2:19], list(max_df$recombination, max_df$hard_sweep), sd) sd_min_df <- aggregate(min_df[, 2:19], list(min_df$recombination, min_df$hard_sweep), sd) colnames(mean_max_df)[1:2] <- c("recombination", "hard_sweep") colnames(mean_min_df)[1:2] <- c("recombination", "hard_sweep") colnames(sd_max_df)[1:2] <- c("recombination", "hard_sweep") colnames(sd_min_df)[1:2] <- c("recombination", "hard_sweep") mean_acc_df <- mean_max_df sd_acc_df <- sd_max_df mean_acc_df <- data.frame(mean_acc_df, "recombination_sweep" = paste(mean_acc_df$recombination, mean_acc_df$hard_sweep, sep = "_")) sd_acc_df <- data.frame(sd_acc_df, "recombination_sweep" = paste(sd_acc_df$recombination, sd_acc_df$hard_sweep, sep = "_")) mean_acc_df$recombination_sweep <- str_replace(mean_acc_df$recombination_sweep, "FALSE", "soft") mean_acc_df$recombination_sweep <- str_replace(mean_acc_df$recombination_sweep, "TRUE", "hard") sd_acc_df$recombination_sweep <- str_replace(sd_acc_df$recombination_sweep, "FALSE", "soft") sd_acc_df$recombination_sweep <- str_replace(sd_acc_df$recombination_sweep, "TRUE", "hard") # determining if lower or higher values have higher predictive accuracy for (i in 1:nrow(mean_max_df)) { for (j in 3:ncol(mean_max_df)) { if (mean_min_df[i, j] > mean_max_df[i, j]) { mean_acc_df[i, j] <- mean_min_df[i, j] sd_acc_df[i, j] <- sd_min_df[i, j] } } } all_acc_df <- max_df all_acc_df <- data.frame(all_acc_df, "recombination_sweep" = paste(all_acc_df$recombination, all_acc_df$hard_sweep, sep = "_")) all_acc_df$recombination_sweep <- str_replace(all_acc_df$recombination_sweep, "FALSE", "soft") all_acc_df$recombination_sweep <- str_replace(all_acc_df$recombination_sweep, "TRUE", "hard") for (i in 1:nrow(max_df)) { for (j in 2:(ncol(max_df) - 2)) { if (min_df[i, j] > max_df[i, j]) { all_acc_df[i, j] <- min_df[i, j] } } } mean_acc_df <- mean_acc_df[, -c(1, 2)] sd_acc_df <- sd_acc_df[, -c(1, 2)] write.table(mean_acc_df, "rev_mean_accuracy_18_prediction_1in1_v2.tsv", quote = FALSE, sep = "\t", row.names = FALSE) write.table(sd_acc_df, "rev_sd_accuracy_18_prediction_1in1_v2.tsv", quote = FALSE, sep = "\t", row.names = FALSE) write.table(all_acc_df, "rev_accuracy_18_prediction_1in1_v2.tsv", quote = FALSE, sep = "\t", row.names = FALSE) ### ### neutral simulations & significance calculations ### ### glob_lengths <- rep(NA, length(title_vec_neutral)) # initialising data frames maxp1_freq_df_neutral <- as.data.frame(matrix(rep(0, 18 * length(title_vec_neutral)), length(title_vec_neutral), 18)) maxp2_freq_df_neutral <- as.data.frame(matrix(rep(0, 18 * length(title_vec_neutral)), length(title_vec_neutral), 18)) minp1_freq_df_neutral <- as.data.frame(matrix(rep(0, 18 * length(title_vec_neutral)), length(title_vec_neutral), 18)) minp2_freq_df_neutral <- as.data.frame(matrix(rep(0, 18 * length(title_vec_neutral)), length(title_vec_neutral), 18)) maxp1_df_neutral <- as.data.frame(matrix(rep(0, 18 * length(title_vec_neutral)), length(title_vec_neutral), 18)) maxp2_df_neutral <- as.data.frame(matrix(rep(0, 18 * length(title_vec_neutral)), length(title_vec_neutral), 18)) minp1_df_neutral <- as.data.frame(matrix(rep(0, 18 * length(title_vec_neutral)), length(title_vec_neutral), 18)) minp2_df_neutral <- as.data.frame(matrix(rep(0, 18 * length(title_vec_neutral)), length(title_vec_neutral), 18)) sig_p_frame <- as.data.frame(matrix(rep(0, 18 * length(title_vec_neutral) * 4), length(title_vec_neutral) * 4, 18)) sig_FDR_frame <- as.data.frame(matrix(rep(0, 18 * length(title_vec_neutral) * 4), length(title_vec_neutral) * 4, 18)) NA_frame <- as.data.frame(matrix(rep(0, 18 * length(title_vec_neutral) * 4), length(title_vec_neutral) * 4, 18)) # looping through all scenarios for (i in 1:length(title_vec_neutral)) { print(i) glob_vec_neutral <- Sys.glob(file.path(paste(path3, scenario_vec_neutral[i], sep = ""))) glob_vec_p3_neutral <- Sys.glob(file.path(paste(path3, scenario_vec_p3_neutral[i], sep = ""))) glob_vec2_neutral <- Sys.glob(file.path(paste(path3, fst_vec_neutral[i], sep = ""))) glob_vec3_neutral <- Sys.glob(file.path(paste(path4, beta_vec_neutral[i], sep = ""))) glob_vec4_neutral <- Sys.glob(file.path(paste(path5, hscan_vec_neutral[i], ".windows", sep = ""))) fv1_frame_neutral <- data.frame(matrix(nrow = 0, ncol = 19)) fv2_frame_neutral <- data.frame(matrix(nrow = 0, ncol = 19)) for (j in 1:100) { # reading in all metric input files fv_neutral <- read.table(glob_vec_neutral[j], header = TRUE) fv_p3_neutral <- read.table(glob_vec_p3_neutral[j], header = TRUE) statvars_neutral <- colsplit(colnames(fv_neutral), "_win", c("metric", "window")) fv_neutral <- unlist(sapply(fv_neutral, unlist)) fv_p3_neutral <- unlist(sapply(fv_p3_neutral, unlist)) fv_neutral <- as.data.frame(matrix(fv_neutral, nrow = 50)) fv_p3_neutral <- as.data.frame(matrix(fv_p3_neutral, nrow = 50)) rownames(fv_neutral) <- 1:50 rownames(fv_p3_neutral) <- 1:50 colnames(fv_neutral) <- unique(statvars_neutral[, 1]) colnames(fv_p3_neutral) <- unique(statvars_neutral[, 1]) fst_frame_neutral <- read.table(glob_vec2_neutral[j], header = TRUE) hscan_frame_neutral <- read.table(glob_vec4_neutral[j], header = TRUE) fv_neutral <- data.frame(fv_neutral, "hscan" = hscan_frame_neutral[1:50, 2]) fv_neutral <- data.frame(fv_neutral, "delta_pi" = fv_neutral[1:50, 1] - fv_p3_neutral[1:50, 1]) fv_neutral <- data.frame(fv_neutral, "Fst" = fst_frame_neutral[1:50, 6]) beta_frame_neutral <- readRDS(glob_vec3_neutral[j]) rownames(beta_frame_neutral) <- 1:50 fv_neutral <- data.frame(fv_neutral, "beta_q0" = beta_frame_neutral[1:50, 2]) fv_neutral <- data.frame(fv_neutral, "beta_q1" = beta_frame_neutral[1:50, 3]) fv_neutral <- data.frame(fv_neutral, "beta_q2" = beta_frame_neutral[1:50, 4]) fv_neutral <- data.frame(fv_neutral, "window_num" = (1:50) * 4000050 / 50) fv1_neutral <- fv_neutral[1:25, ] fv2_neutral <- fv_neutral[26:50, ] colnames(fv1_frame_neutral) <- colnames(fv1_neutral) colnames(fv2_frame_neutral) <- colnames(fv2_neutral) for (k in 1:(ncol(fv_neutral) - 1)) { part1_neutral <- fv1_neutral[, k] part2_neutral <- fv2_neutral[, k] maxp1_1_neutral <- order(part1_neutral, decreasing = TRUE)[1] maxp1_2_neutral <- order(part1_neutral, decreasing = TRUE)[2] maxp2_1_neutral <- order(part2_neutral, decreasing = TRUE)[1] maxp2_2_neutral <- order(part2_neutral, decreasing = TRUE)[2] minp1_1_neutral <- order(part1_neutral, decreasing = FALSE)[1] minp1_2_neutral <- order(part1_neutral, decreasing = FALSE)[2] minp2_1_neutral <- order(part2_neutral, decreasing = FALSE)[1] minp2_2_neutral <- order(part2_neutral, decreasing = FALSE)[2] #if ((maxp1_1_neutral %in% 20:21) && (maxp1_2_neutral %in% 20:21)) { if (maxp1_1_neutral == 13) { maxp1_df_neutral[i, k] <- maxp1_df_neutral[i, k] + 1 } #if ((maxp2_1_neutral %in% 60:61) && (maxp2_2_neutral %in% 60:61)) { if (maxp2_1_neutral == 13) { maxp2_df_neutral[i, k] <- maxp2_df_neutral[i, k] + 1 } #if ((minp1_1 %in% 20:21) && (minp1_2 %in% 20:21)) { if (minp1_1_neutral == 13) { minp1_df_neutral[i, k] <- minp1_df_neutral[i, k] + 1 } #if ((minp2_1 %in% 60:61) && (minp2_2 %in% 60:61)) { if (minp2_1_neutral == 13) { minp2_df_neutral[i, k] <- minp2_df_neutral[i, k] + 1 } } fv1_frame_neutral <- rbind(fv1_frame_neutral, fv1_neutral[13, ]) fv2_frame_neutral <- rbind(fv2_frame_neutral, fv2_neutral[13, ]) } ia <- i + (floor((i + 3) / 4) * 4) - 4 ib <- i + (floor((i + 3) / 4) * 4) glob_vec_a <- Sys.glob(file.path(paste(scenario_vec[ia], sep = ""))) glob_vec_a_p3 <- Sys.glob(file.path(paste(scenario_vec_p3[ia], sep = ""))) glob_vec_b <- Sys.glob(file.path(paste(scenario_vec[ib], sep = ""))) glob_vec_b_p3 <- Sys.glob(file.path(paste(scenario_vec_p3[ib], sep = ""))) glob_vec2_a <- Sys.glob(file.path(paste(path, fst_vec[ia], sep = ""))) glob_vec2_b <- Sys.glob(file.path(paste(path, fst_vec[ib], sep = ""))) glob_vec3_a <- Sys.glob(file.path(paste(path2, beta_vec[ia], sep = ""))) glob_vec3_b <- Sys.glob(file.path(paste(path2, beta_vec[ib], sep = ""))) glob_vec4_a <- Sys.glob(file.path(paste(path, hscan_vec[ia], ".windows", sep = ""))) glob_vec4_b <- Sys.glob(file.path(paste(path, hscan_vec[ib], ".windows", sep = ""))) pval_list <- list() pval_list2 <- list() # 1: glob_vec_a for (j in 1:100) { ### print("g") curr_frame_a <- read.table(glob_vec_a[j], header = TRUE) print("g2") curr_frame_a <- unlist(curr_frame_a) curr_frame_a <- as.data.frame(matrix(curr_frame_a, nrow = 50)) curr_frame_a2 <- curr_frame_a[38, ] curr_frame_a <- curr_frame_a[13, ] ### for (k in 1:12) { iscen <- 1 if (mean_min_df[iscen, k + 2] < mean_max_df[iscen, k + 2]) { if ((!is.na(curr_frame_a[, k])) && (!is.nan(curr_frame_a[, k]))) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval <- t.test(fv1_frame_neutral[, k], alternative = "less", mu = curr_frame_a[, k])$p.value } else { curr_pval <- NA } } else { if ((!is.na(curr_frame_a[, k])) && (!is.nan(curr_frame_a[, k]))) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval <- t.test(fv1_frame_neutral[, k], alternative = "greater", mu = curr_frame_a[, k])$p.value } else { curr_pval <- NA } } if (mean_min_df[iscen + 1, k + 2] < mean_max_df[iscen + 1, k + 2]) { if ((!is.na(curr_frame_a2[, k])) && (!is.nan(curr_frame_a2[, k]))) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval2 <- t.test(fv2_frame_neutral[, k], alternative = "less", mu = curr_frame_a2[, k])$p.value } else { curr_pval2 <- NA } } else { if ((!is.na(curr_frame_a2[, k])) && (!is.nan(curr_frame_a2[, k]))) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval2 <- t.test(fv2_frame_neutral[, k], alternative = "greater", mu = curr_frame_a2[, k])$p.value } else { curr_pval2 <- NA } } if (length(pval_list) >= k) { pval_list[[k]] <- c(pval_list[[k]], curr_pval) } else { pval_list[[k]] <- curr_pval } if (length(pval_list2) >= k) { pval_list2[[k]] <- c(pval_list2[[k]], curr_pval2) } else { pval_list2[[k]] <- curr_pval2 } if (!is.na(curr_pval)) { if (curr_pval < 0.05) { sig_p_frame[ia, k] <- sig_p_frame[ia, k] + 1 } } else { NA_frame[ia, k] <- NA_frame[ia, k] + 1 } if (!is.na(curr_pval2)) { if (curr_pval2 < 0.05) { sig_p_frame[ia + 96, k] <- sig_p_frame[ia + 96, k] + 1 } } else { NA_frame[ia + 96, k] <- NA_frame[ia + 96, k] + 1 } } } # 1: FDRs FDR_list <- lapply(pval_list, function(x) p.adjust(x, method = "fdr")) FDR_list2 <- lapply(pval_list2, function(x) p.adjust(x, method = "fdr")) for (k in 1:12) { for (l in 1:length(FDR_list[[k]])) { if (!is.na(FDR_list[[k]][l])) { if (FDR_list[[k]][l] < 0.1) { sig_FDR_frame[ia, k] <- sig_FDR_frame[ia, k] + 1 } } } } for (k in 1:12) { for (l in 1:length(FDR_list2[[k]])) { if (!is.na(FDR_list2[[k]][l])) { if (FDR_list2[[k]][l] < 0.1) { sig_FDR_frame[ia + 96, k] <- sig_FDR_frame[ia + 96, k] + 1 } } } } pval_list <- list() pval_list2 <- list() # 1 deltas: glob_vec_a for (j in 1:100) { ### curr_frame_a <- read.table(glob_vec_a[j], header = TRUE) curr_frame_a_p3 <- read.table(glob_vec_a_p3[j], header = TRUE) curr_frame_a <- unlist(sapply(curr_frame_a, unlist)) curr_frame_a_p3 <- unlist(sapply(curr_frame_a_p3, unlist)) curr_frame_a <- as.data.frame(matrix(curr_frame_a, nrow = 50)) curr_frame_a_p3 <- as.data.frame(matrix(curr_frame_a_p3, nrow = 50)) curr_frame_a2 <- curr_frame_a[38, ] - curr_frame_a_p3[38, ] curr_frame_a <- curr_frame_a[13, ] - curr_frame_a_p3[13, ] ### for (k in 13) { iscen <- 1 if (mean_min_df[iscen, k + 2] < mean_max_df[iscen, k + 2]) { if ((!is.na(curr_frame_a[, k - 12])) && (!is.nan(curr_frame_a[, k - 12]))) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval <- t.test(fv1_frame_neutral[, k + 1], alternative = "less", mu = curr_frame_a[, k - 12])$p.value } else { curr_pval <- NA } } else { if ((!is.na(curr_frame_a[, k - 12])) && (!is.nan(curr_frame_a[, k - 12]))) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval <- t.test(fv1_frame_neutral[, k + 1], alternative = "greater", mu = curr_frame_a[, k - 12])$p.value } else { curr_pval <- NA } } if (mean_min_df[iscen + 1, k + 2] < mean_max_df[iscen + 1, k + 2]) { if ((!is.na(curr_frame_a2[, k - 12])) && (!is.nan(curr_frame_a2[, k - 12]))) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval2 <- t.test(fv2_frame_neutral[, k + 1], alternative = "less", mu = curr_frame_a2[, k - 12])$p.value } else { curr_pval2 <- NA } } else { if ((!is.na(curr_frame_a2[, k - 12])) && (!is.nan(curr_frame_a2[, k - 12]))) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval2 <- t.test(fv2_frame_neutral[, k + 1], alternative = "greater", mu = curr_frame_a2[, k - 12])$p.value } else { curr_pval2 <- NA } } if (length(pval_list) >= k) { pval_list[[k]] <- c(pval_list[[k]], curr_pval) } else { pval_list[[k]] <- curr_pval } if (length(pval_list2) >= k) { pval_list2[[k]] <- c(pval_list2[[k]], curr_pval2) } else { pval_list2[[k]] <- curr_pval2 } if (!is.na(curr_pval)) { if (curr_pval < 0.05) { sig_p_frame[ia, k] <- sig_p_frame[ia, k] + 1 } } else { NA_frame[ia, k] <- NA_frame[ia, k] + 1 } if (!is.na(curr_pval2)) { if (curr_pval2 < 0.05) { sig_p_frame[ia + 96, k] <- sig_p_frame[ia + 96, k] + 1 } } else { NA_frame[ia + 96, k] <- NA_frame[ia + 96, k] + 1 } } } # 1 deltas: FDRs FDR_list <- lapply(pval_list, function(x) p.adjust(x, method = "fdr")) FDR_list2 <- lapply(pval_list2, function(x) p.adjust(x, method = "fdr")) for (k in 13) { for (l in 1:length(FDR_list[[k]])) { if (!is.na(FDR_list[[k]][l])) { if (FDR_list[[k]][l] < 0.1) { sig_FDR_frame[ia, k] <- sig_FDR_frame[ia, k] + 1 } } } } for (k in 13) { for (l in 1:length(FDR_list2[[k]])) { if (!is.na(FDR_list2[[k]][l])) { if (FDR_list2[[k]][l] < 0.1) { sig_FDR_frame[ia + 96, k] <- sig_FDR_frame[ia + 96, k] + 1 } } } } pval_list <- list() pval_list2 <- list() # 2: glob_vec_b for (j in 1:100) { ### curr_frame_a <- read.table(glob_vec_b[j], header = TRUE) curr_frame_a <- unlist(sapply(curr_frame_a, unlist)) curr_frame_a <- as.data.frame(matrix(curr_frame_a, nrow = 50)) curr_frame_a2 <- curr_frame_a[38, ] curr_frame_a <- curr_frame_a[13, ] ### for (k in 1:12) { iscen <- 3 if (mean_min_df[iscen, k + 2] < mean_max_df[iscen, k + 2]) { if ((!is.na(curr_frame_a[, k])) && (!is.nan(curr_frame_a[, k]))) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval <- t.test(fv1_frame_neutral[, k], alternative = "less", mu = curr_frame_a[, k])$p.value } else { curr_pval <- NA } } else { if ((!is.na(curr_frame_a[, k])) && (!is.nan(curr_frame_a[, k]))) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval <- t.test(fv1_frame_neutral[, k], alternative = "greater", mu = curr_frame_a[, k])$p.value } else { curr_pval <- NA } } if (mean_min_df[iscen + 1, k + 2] < mean_max_df[iscen + 1, k + 2]) { if ((!is.na(curr_frame_a2[, k])) && (!is.nan(curr_frame_a2[, k]))) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval2 <- t.test(fv2_frame_neutral[, k], alternative = "less", mu = curr_frame_a2[, k])$p.value } else { curr_pval2 <- NA } } else { if ((!is.na(curr_frame_a2[, k])) && (!is.nan(curr_frame_a2[, k]))) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval2 <- t.test(fv2_frame_neutral[, k], alternative = "greater", mu = curr_frame_a2[, k])$p.value } else { curr_pval2 <- NA } } if (length(pval_list) >= k) { pval_list[[k]] <- c(pval_list[[k]], curr_pval) } else { pval_list[[k]] <- curr_pval } if (length(pval_list2) >= k) { pval_list2[[k]] <- c(pval_list2[[k]], curr_pval2) } else { pval_list2[[k]] <- curr_pval2 } if (!is.na(curr_pval)) { if (curr_pval < 0.05) { sig_p_frame[ib, k] <- sig_p_frame[ib, k] + 1 } } else { NA_frame[ib, k] <- NA_frame[ib, k] + 1 } if (!is.na(curr_pval2)) { if (curr_pval2 < 0.05) { sig_p_frame[ib + 96, k] <- sig_p_frame[ib + 96, k] + 1 } } else { NA_frame[ib + 96, k] <- NA_frame[ib + 96, k] + 1 } } } # 2: FDRs FDR_list <- lapply(pval_list, function(x) p.adjust(x, method = "fdr")) FDR_list2 <- lapply(pval_list2, function(x) p.adjust(x, method = "fdr")) for (k in 1:12) { for (l in 1:length(FDR_list[[k]])) { if (!is.na(FDR_list[[k]][l])) { if (FDR_list[[k]][l] < 0.1) { sig_FDR_frame[ib, k] <- sig_FDR_frame[ib, k] + 1 } } } } for (k in 1:12) { for (l in 1:length(FDR_list2[[k]])) { if (!is.na(FDR_list2[[k]][l])) { if (FDR_list2[[k]][l] < 0.1) { sig_FDR_frame[ib + 96, k] <- sig_FDR_frame[ib + 96, k] + 1 } } } } pval_list <- list() pval_list2 <- list() # 2 deltas: glob_vec_b for (j in 1:100) { ### curr_frame_a <- read.table(glob_vec_b[j], header = TRUE) curr_frame_a_p3 <- read.table(glob_vec_b_p3[j], header = TRUE) curr_frame_a <- unlist(sapply(curr_frame_a, unlist)) curr_frame_a_p3 <- unlist(sapply(curr_frame_a_p3, unlist)) curr_frame_a <- as.data.frame(matrix(curr_frame_a, nrow = 50)) curr_frame_a_p3 <- as.data.frame(matrix(curr_frame_a_p3, nrow = 50)) curr_frame_a2 <- curr_frame_a[38, ] - curr_frame_a_p3[38, ] curr_frame_a <- curr_frame_a[13, ] - curr_frame_a_p3[13, ] ### for (k in 13) { iscen <- 3 if (mean_min_df[iscen, k + 2] < mean_max_df[iscen, k + 2]) { if ((!is.na(curr_frame_a[, k - 12])) && (!is.nan(curr_frame_a[, k - 12]))) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval <- t.test(fv1_frame_neutral[, k + 1], alternative = "less", mu = curr_frame_a[, k - 12])$p.value } else { curr_pval <- NA } } else { if ((!is.na(curr_frame_a[, k - 12])) && (!is.nan(curr_frame_a[, k - 12]))) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval <- t.test(fv1_frame_neutral[, k + 1], alternative = "greater", mu = curr_frame_a[, k - 12])$p.value } else { curr_pval <- NA } } if (mean_min_df[iscen + 1, k + 2] < mean_max_df[iscen + 1, k + 2]) { if ((!is.na(curr_frame_a2[, k - 12])) && (!is.nan(curr_frame_a2[, k - 12]))) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval2 <- t.test(fv2_frame_neutral[, k + 1], alternative = "less", mu = curr_frame_a2[, k - 12])$p.value } else { curr_pval2 <- NA } } else { if ((!is.na(curr_frame_a2[, k - 12])) && (!is.nan(curr_frame_a2[, k - 12]))) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval2 <- t.test(fv2_frame_neutral[, k + 1], alternative = "greater", mu = curr_frame_a2[, k - 12])$p.value } else { curr_pval2 <- NA } } if (length(pval_list) >= k) { pval_list[[k]] <- c(pval_list[[k]], curr_pval) } else { pval_list[[k]] <- curr_pval } if (length(pval_list2) >= k) { pval_list2[[k]] <- c(pval_list2[[k]], curr_pval2) } else { pval_list2[[k]] <- curr_pval2 } if (!is.na(curr_pval)) { if (curr_pval < 0.05) { sig_p_frame[ib, k] <- sig_p_frame[ib, k] + 1 } } else { NA_frame[ib, k] <- NA_frame[ib, k] + 1 } if (!is.na(curr_pval2)) { if (curr_pval2 < 0.05) { sig_p_frame[ib + 96, k] <- sig_p_frame[ib + 96, k] + 1 } } else { NA_frame[ib + 96, k] <- NA_frame[ib + 96, k] + 1 } } } # 2 deltas: FDRs FDR_list <- lapply(pval_list, function(x) p.adjust(x, method = "fdr")) FDR_list2 <- lapply(pval_list2, function(x) p.adjust(x, method = "fdr")) for (k in 13) { for (l in 1:length(FDR_list[[k]])) { if (!is.na(FDR_list[[k]][l])) { if (FDR_list[[k]][l] < 0.1) { sig_FDR_frame[ib, k] <- sig_FDR_frame[ib, k] + 1 } } } } for (k in 13) { for (l in 1:length(FDR_list2[[k]])) { if (!is.na(FDR_list2[[k]][l])) { if (FDR_list2[[k]][l] < 0.1) { sig_FDR_frame[ib + 96, k] <- sig_FDR_frame[ib + 96, k] + 1 } } } } pval_list <- list() pval_list2 <- list() # 3: glob_vec2_a for (j in 1:100) { ### curr_frame_a <- read.table(glob_vec2_a[j], header = TRUE) curr_frame_a <- unlist(sapply(curr_frame_a, unlist)) curr_frame_a <- as.data.frame(matrix(curr_frame_a, nrow = 50)) if (any(!seq(1, 3920050, 80001) %in% curr_frame_a[, 2])) { for (j2 in 1:length(which(!seq(1, 3920050, 80001) %in% curr_frame_a[, 2]))) { missing_row <- seq(1, 3920050, 80001)[j2] if (missing_row %in% seq(80002, 3840049, 80001)) { position <- which(missing_row %in% seq(1, 3920050, 80001)) new_row <- c(1, missing_row, missing_row + 80000, 0, mean(curr_frame_a[position - 1, 5], curr_frame_a[position + 1, 5]), mean(curr_frame_a[position - 1, 6], curr_frame_a[position + 1, 6])) curr_frame_a <- rbind(curr_frame_a[1:(position - 1), ], new_row, curr_frame_a[(position + 1):nrow(curr_frame_a), ]) } else if (missing_row == 1) { position <- 1 new_row <- c(1, missing_row, missing_row + 80000, 0, curr_frame_a[position + 1, 5], curr_frame_a[position + 1, 6]) curr_frame_a <- rbind(new_row, curr_frame_a[(position + 1):nrow(curr_frame_a), ]) } else { position <- 50 new_row <- c(1, missing_row, missing_row + 80000, 0, curr_frame_a[position - 1, 5], curr_frame_a[position - 1, 6]) curr_frame_a <- rbind(curr_frame_a[(position - 1):nrow(curr_frame_a), ], new_row) } } } curr_frame_a <- curr_frame_a[1:50, ] curr_frame_a2 <- curr_frame_a[38, ] curr_frame_a <- curr_frame_a[13, ] ### for (k in 14) { iscen <- 1 if (mean_min_df[iscen, k + 2] < mean_max_df[iscen, k + 2]) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) if ((!is.na(curr_frame_a[, ncol(curr_frame_a)])) && (!is.nan(curr_frame_a[, ncol(curr_frame_a)])) && (ncol(curr_frame_a) == 6)) { curr_pval <- t.test(fv1_frame_neutral[, k + 1], alternative = "less", mu = curr_frame_a[, 6])$p.value } else { curr_pval <- NA } } else { if ((!is.na(curr_frame_a[, ncol(curr_frame_a)])) && (!is.nan(curr_frame_a[, ncol(curr_frame_a)])) && (ncol(curr_frame_a) == 6)) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval <- t.test(fv1_frame_neutral[, k + 1], alternative = "greater", mu = curr_frame_a[, 6])$p.value } else { curr_pval <- NA } } if (mean_min_df[iscen + 1, k + 2] < mean_max_df[iscen + 1, k + 2]) { if ((!is.na(curr_frame_a2[, ncol(curr_frame_a2)])) && (!is.nan(curr_frame_a2[, ncol(curr_frame_a2)])) && (ncol(curr_frame_a2) == 6)) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval2 <- t.test(fv2_frame_neutral[, k + 1], alternative = "less", mu = curr_frame_a2[, 6])$p.value } else { curr_pval2 <- NA } } else { if ((!is.na(curr_frame_a2[, ncol(curr_frame_a2)])) && (!is.nan(curr_frame_a2[, ncol(curr_frame_a2)])) && (ncol(curr_frame_a2) == 6)) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval2 <- t.test(fv2_frame_neutral[, k + 1], alternative = "greater", mu = curr_frame_a2[, 6])$p.value } else { curr_pval2 <- NA } } if (length(pval_list) >= k) { pval_list[[k]] <- c(pval_list[[k]], curr_pval) } else { pval_list[[k]] <- curr_pval } if (length(pval_list2) >= k) { pval_list2[[k]] <- c(pval_list2[[k]], curr_pval2) } else { pval_list2[[k]] <- curr_pval2 } if (!is.na(curr_pval)) { if (curr_pval < 0.05) { sig_p_frame[ia, k] <- sig_p_frame[ia, k] + 1 } } else { NA_frame[ia, k] <- NA_frame[ia, k] + 1 } if (!is.na(curr_pval2)) { if (curr_pval2 < 0.05) { sig_p_frame[ia + 96, k] <- sig_p_frame[ia + 96, k] + 1 } } else { NA_frame[ia + 96, k] <- NA_frame[ia + 96, k] + 1 } } } # 3: FDRs FDR_list <- lapply(pval_list, function(x) p.adjust(x, method = "fdr")) FDR_list2 <- lapply(pval_list2, function(x) p.adjust(x, method = "fdr")) for (k in 14) { for (l in 1:length(FDR_list[[k]])) { if (!is.na(FDR_list[[k]][l])) { if (FDR_list[[k]][l] < 0.1) { sig_FDR_frame[ia, k] <- sig_FDR_frame[ia, k] + 1 } } } } for (k in 14) { for (l in 1:length(FDR_list2[[k]])) { if (!is.na(FDR_list2[[k]][l])) { if (FDR_list2[[k]][l] < 0.1) { sig_FDR_frame[ia + 96, k] <- sig_FDR_frame[ia + 96, k] + 1 } } } } pval_list <- list() pval_list2 <- list() # 4: glob_vec2_b for (j in 1:100) { ### curr_frame_a <- read.table(glob_vec2_b[j], header = TRUE) curr_frame_a <- unlist(sapply(curr_frame_a, unlist)) curr_frame_a <- as.data.frame(matrix(curr_frame_a, nrow = 50)) if (any(!seq(1, 3920050, 80001) %in% curr_frame_a[, 2])) { for (j2 in 1:length(which(!seq(1, 3920050, 80001) %in% curr_frame_a[, 2]))) { missing_row <- seq(1, 3920050, 80001)[j2] if (missing_row %in% seq(80002, 3840049, 80001)) { position <- which(missing_row %in% seq(1, 3920050, 80001)) new_row <- c(1, missing_row, missing_row + 80000, 0, mean(curr_frame_a[position - 1, 5], curr_frame_a[position + 1, 5]), mean(curr_frame_a[position - 1, 6], curr_frame_a[position + 1, 6])) curr_frame_a <- rbind(curr_frame_a[1:(position - 1), ], new_row, curr_frame_a[(position + 1):nrow(curr_frame_a), ]) } else if (missing_row == 1) { position <- 1 new_row <- c(1, missing_row, missing_row + 80000, 0, curr_frame_a[position + 1, 5], curr_frame_a[position + 1, 6]) curr_frame_a <- rbind(new_row, curr_frame_a[(position + 1):nrow(curr_frame_a), ]) } else { position <- 50 new_row <- c(1, missing_row, missing_row + 80000, 0, curr_frame_a[position - 1, 5], curr_frame_a[position - 1, 6]) curr_frame_a <- rbind(curr_frame_a[(position - 1):nrow(curr_frame_a), ], new_row) } } } curr_frame_a <- curr_frame_a[1:50, ] curr_frame_a2 <- curr_frame_a[38, ] curr_frame_a <- curr_frame_a[13, ] ### for (k in 14) { iscen <- 3 if (mean_min_df[iscen, k + 2] < mean_max_df[iscen, k + 2]) { if ((!is.na(curr_frame_a[, ncol(curr_frame_a)])) && (!is.nan(curr_frame_a[, ncol(curr_frame_a)])) && (ncol(curr_frame_a) == 6)) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval <- t.test(fv1_frame_neutral[, k + 1], alternative = "less", mu = curr_frame_a[, 6])$p.value } else { curr_pval <- NA } } else { if ((!is.na(curr_frame_a[, ncol(curr_frame_a)])) && (!is.nan(curr_frame_a[, ncol(curr_frame_a)])) && (ncol(curr_frame_a) == 6)) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval <- t.test(fv1_frame_neutral[, k + 1], alternative = "greater", mu = curr_frame_a[, 6])$p.value } else { curr_pval <- NA } } if (mean_min_df[iscen + 1, k + 2] < mean_max_df[iscen + 1, k + 2]) { if ((!is.na(curr_frame_a2[, ncol(curr_frame_a2)])) && (!is.nan(curr_frame_a2[, ncol(curr_frame_a2)])) && (ncol(curr_frame_a2) == 6)) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval2 <- t.test(fv2_frame_neutral[, k + 1], alternative = "less", mu = curr_frame_a2[, 6])$p.value } else { curr_pval2 <- NA } } else { if ((!is.na(curr_frame_a2[, ncol(curr_frame_a2)])) && (!is.nan(curr_frame_a2[, ncol(curr_frame_a2)])) && (ncol(curr_frame_a2) == 6)) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval2 <- t.test(fv2_frame_neutral[, k + 1], alternative = "greater", mu = curr_frame_a2[, 6])$p.value } else { curr_pval2 <- NA } } if (length(pval_list) >= k) { pval_list[[k]] <- c(pval_list[[k]], curr_pval) } else { pval_list[[k]] <- curr_pval } if (length(pval_list2) >= k) { pval_list2[[k]] <- c(pval_list2[[k]], curr_pval2) } else { pval_list2[[k]] <- curr_pval2 } if (!is.na(curr_pval)) { if (curr_pval < 0.05) { sig_p_frame[ib, k] <- sig_p_frame[ib, k] + 1 } } else { NA_frame[ib, k] <- NA_frame[ib, k] + 1 } if (!is.na(curr_pval2)) { if (curr_pval2 < 0.05) { sig_p_frame[ib + 96, k] <- sig_p_frame[ib + 96, k] + 1 } } else { NA_frame[ib + 96, k] <- NA_frame[ib + 96, k] + 1 } } } # 4: FDRs FDR_list <- lapply(pval_list, function(x) p.adjust(x, method = "fdr")) FDR_list2 <- lapply(pval_list2, function(x) p.adjust(x, method = "fdr")) for (k in 14) { for (l in 1:length(FDR_list[[k]])) { if (!is.na(FDR_list[[k]][l])) { if (FDR_list[[k]][l] < 0.1) { sig_FDR_frame[ib, k] <- sig_FDR_frame[ib, k] + 1 } } } } for (k in 14) { for (l in 1:length(FDR_list2[[k]])) { if (!is.na(FDR_list2[[k]][l])) { if (FDR_list2[[k]][l] < 0.1) { sig_FDR_frame[ib + 96, k] <- sig_FDR_frame[ib + 96, k] + 1 } } } } pval_list <- list() pval_list2 <- list() # 5: glob_vec4_a for (j in 1:100) { ### curr_frame_a <- read.table(glob_vec4_a[j], header = TRUE) curr_frame_a <- unlist(sapply(curr_frame_a, unlist)) curr_frame_a <- as.data.frame(matrix(curr_frame_a, nrow = 50)) curr_frame_a <- curr_frame_a[1:50, ] curr_frame_a2 <- curr_frame_a[38, ] curr_frame_a <- curr_frame_a[13, ] ### for (k in 15) { iscen <- 1 if (mean_min_df[iscen, k + 2] < mean_max_df[iscen, k + 2]) { if ((!is.na(curr_frame_a[, ncol(curr_frame_a)])) && (!is.nan(curr_frame_a[, ncol(curr_frame_a)])) && (ncol(curr_frame_a) == 2)) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval <- t.test(fv1_frame_neutral[, k - 2], alternative = "less", mu = curr_frame_a[, 2])$p.value } else { curr_pval <- NA } } else { if ((!is.na(curr_frame_a[, ncol(curr_frame_a)])) && (!is.nan(curr_frame_a[, ncol(curr_frame_a)])) && (ncol(curr_frame_a) == 2)) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval <- t.test(fv1_frame_neutral[, k - 2], alternative = "greater", mu = curr_frame_a[, 2])$p.value } else { curr_pval <- NA } } if (mean_min_df[iscen + 1, k + 2] < mean_max_df[iscen + 1, k + 2]) { if ((!is.na(curr_frame_a2[, ncol(curr_frame_a2)])) && (!is.nan(curr_frame_a2[, ncol(curr_frame_a2)])) && (ncol(curr_frame_a2) == 2)) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval2 <- t.test(fv2_frame_neutral[, k - 2], alternative = "less", mu = curr_frame_a2[, 2])$p.value } else { curr_pval2 <- NA } } else { if ((!is.na(curr_frame_a2[, ncol(curr_frame_a2)])) && (!is.nan(curr_frame_a2[, ncol(curr_frame_a2)])) && (ncol(curr_frame_a2) == 2)) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval2 <- t.test(fv2_frame_neutral[, k - 2], alternative = "greater", mu = curr_frame_a2[, 2])$p.value } else { curr_pval2 <- NA } } if (length(pval_list) >= k) { pval_list[[k]] <- c(pval_list[[k]], curr_pval) } else { pval_list[[k]] <- curr_pval } if (length(pval_list2) >= k) { pval_list2[[k]] <- c(pval_list2[[k]], curr_pval2) } else { pval_list2[[k]] <- curr_pval2 } if (!is.na(curr_pval)) { if (curr_pval < 0.05) { sig_p_frame[ia, k] <- sig_p_frame[ia, k] + 1 } } else { NA_frame[ia, k] <- NA_frame[ia, k] + 1 } if (!is.na(curr_pval2)) { if (curr_pval2 < 0.05) { sig_p_frame[ia + 96, k] <- sig_p_frame[ia + 96, k] + 1 } } else { NA_frame[ia + 96, k] <- NA_frame[ia + 96, k] + 1 } } } # 5: FDRs FDR_list <- lapply(pval_list, function(x) p.adjust(x, method = "fdr")) FDR_list2 <- lapply(pval_list2, function(x) p.adjust(x, method = "fdr")) for (k in 15) { for (l in 1:length(FDR_list[[k]])) { if (!is.na(FDR_list[[k]][l])) { if (FDR_list[[k]][l] < 0.1) { sig_FDR_frame[ia, k] <- sig_FDR_frame[ia, k] + 1 } } } } for (k in 15) { for (l in 1:length(FDR_list2[[k]])) { if (!is.na(FDR_list2[[k]][l])) { if (FDR_list2[[k]][l] < 0.1) { sig_FDR_frame[ia + 96, k] <- sig_FDR_frame[ia + 96, k] + 1 } } } } pval_list <- list() pval_list2 <- list() # 6: glob_vec4_b for (j in 1:100) { ### curr_frame_a <- read.table(glob_vec4_b[j], header = TRUE) curr_frame_a <- unlist(sapply(curr_frame_a, unlist)) curr_frame_a <- as.data.frame(matrix(curr_frame_a, nrow = 50)) curr_frame_a <- curr_frame_a[1:50, ] curr_frame_a2 <- curr_frame_a[38, ] curr_frame_a <- curr_frame_a[13, ] ### for (k in 15) { iscen <- 1 if (mean_min_df[iscen, k + 2] < mean_max_df[iscen, k + 2]) { if ((!is.na(curr_frame_a[, ncol(curr_frame_a)])) && (!is.nan(curr_frame_a[, ncol(curr_frame_a)])) && (ncol(curr_frame_a) == 2)) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval <- t.test(fv1_frame_neutral[, k - 2], alternative = "less", mu = curr_frame_a[, 2])$p.value } else { curr_pval <- NA } } else { if ((!is.na(curr_frame_a[, ncol(curr_frame_a)])) && (!is.nan(curr_frame_a[, ncol(curr_frame_a)])) && (ncol(curr_frame_a) == 2)) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval <- t.test(fv1_frame_neutral[, k - 2], alternative = "greater", mu = curr_frame_a[, 2])$p.value } else { curr_pval <- NA } } if (mean_min_df[iscen + 1, k + 2] < mean_max_df[iscen + 1, k + 2]) { if ((!is.na(curr_frame_a2[, ncol(curr_frame_a2)])) && (!is.nan(curr_frame_a2[, ncol(curr_frame_a2)])) && (ncol(curr_frame_a2) == 2)) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval2 <- t.test(fv2_frame_neutral[, k - 2], alternative = "less", mu = curr_frame_a2[, 2])$p.value } else { curr_pval2 <- NA } } else { if ((!is.na(curr_frame_a2[, ncol(curr_frame_a2)])) && (!is.nan(curr_frame_a2[, ncol(curr_frame_a2)])) && (ncol(curr_frame_a2) == 2)) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval2 <- t.test(fv2_frame_neutral[, k - 2], alternative = "greater", mu = curr_frame_a2[, 2])$p.value } else { curr_pval2 <- NA } } if (length(pval_list) >= k) { pval_list[[k]] <- c(pval_list[[k]], curr_pval) } else { pval_list[[k]] <- curr_pval } if (length(pval_list2) >= k) { pval_list2[[k]] <- c(pval_list2[[k]], curr_pval2) } else { pval_list2[[k]] <- curr_pval2 } if (!is.na(curr_pval)) { if (curr_pval < 0.05) { sig_p_frame[ib, k] <- sig_p_frame[ib, k] + 1 } } else { NA_frame[ib, k] <- NA_frame[ib, k] + 1 } if (!is.na(curr_pval2)) { if (curr_pval2 < 0.05) { sig_p_frame[ib + 96, k] <- sig_p_frame[ib + 96, k] + 1 } } else { NA_frame[ib + 96, k] <- NA_frame[ib + 96, k] + 1 } } } # 6: FDRs FDR_list <- lapply(pval_list, function(x) p.adjust(x, method = "fdr")) FDR_list2 <- lapply(pval_list2, function(x) p.adjust(x, method = "fdr")) for (k in 15) { for (l in 1:length(FDR_list[[k]])) { if (!is.na(FDR_list[[k]][l])) { if (FDR_list[[k]][l] < 0.1) { sig_FDR_frame[ib, k] <- sig_FDR_frame[ib, k] + 1 } } } } for (k in 15) { for (l in 1:length(FDR_list2[[k]])) { if (!is.na(FDR_list2[[k]][l])) { if (FDR_list2[[k]][l] < 0.1) { sig_FDR_frame[ib + 96, k] <- sig_FDR_frame[ib + 96, k] + 1 } } } } pval_list <- list() pval_list2 <- list() # 7: glob_vec3_a for (j in 1:100) { ### curr_frame_a <- readRDS(glob_vec3_a[j]) curr_frame_a2 <- curr_frame_a[38, 2:4] curr_frame_a <- curr_frame_a[13, 2:4] ### for (k in 16:18) { iscen <- 1 if (mean_min_df[iscen, k + 2] < mean_max_df[iscen, k + 2]) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval <- t.test(fv1_frame_neutral[, k], alternative = "less", mu = curr_frame_a[, k - 15])$p.value } else { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval <- t.test(fv1_frame_neutral[, k], alternative = "greater", mu = curr_frame_a[, k - 15])$p.value } if (mean_min_df[iscen + 1, k + 2] < mean_max_df[iscen + 1, k + 2]) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval2 <- t.test(fv2_frame_neutral[, k], alternative = "less", mu = curr_frame_a2[, k - 15])$p.value } else { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval2 <- t.test(fv2_frame_neutral[, k], alternative = "greater", mu = curr_frame_a2[, k - 15])$p.value } if (length(pval_list) >= k) { pval_list[[k]] <- c(pval_list[[k]], curr_pval) } else { pval_list[[k]] <- curr_pval } if (length(pval_list2) >= k) { pval_list2[[k]] <- c(pval_list2[[k]], curr_pval2) } else { pval_list2[[k]] <- curr_pval2 } if (curr_pval < 0.05) { sig_p_frame[ia, k] <- sig_p_frame[ia, k] + 1 } if (curr_pval2 < 0.05) { sig_p_frame[ia + 96, k] <- sig_p_frame[ia + 96, k] + 1 } } } # 7: FDRs FDR_list <- lapply(pval_list, function(x) p.adjust(x, method = "fdr")) FDR_list2 <- lapply(pval_list2, function(x) p.adjust(x, method = "fdr")) for (k in 16:18) { for (l in 1:length(FDR_list[[k]])) { if (!is.na(FDR_list[[k]][l])) { if (FDR_list[[k]][l] < 0.1) { sig_FDR_frame[ia, k] <- sig_FDR_frame[ia, k] + 1 } } } } for (k in 16:18) { for (l in 1:length(FDR_list2[[k]])) { if (!is.na(FDR_list2[[k]][l])) { if (FDR_list2[[k]][l] < 0.1) { sig_FDR_frame[ia + 96, k] <- sig_FDR_frame[ia + 96, k] + 1 } } } } pval_list <- list() pval_list2 <- list() # 8: glob_vec3_b for (j in 1:100) { ### curr_frame_a <- readRDS(glob_vec3_b[j]) curr_frame_a2 <- curr_frame_a[38, 2:4] curr_frame_a <- curr_frame_a[13, 2:4] ### for (k in 16:18) { iscen <- 3 if (mean_min_df[iscen, k + 2] < mean_max_df[iscen, k + 2]) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval <- t.test(fv1_frame_neutral[, k], alternative = "less", mu = curr_frame_a[, k - 15])$p.value } else { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval <- t.test(fv1_frame_neutral[, k], alternative = "greater", mu = curr_frame_a[, k - 15])$p.value } if (mean_min_df[iscen + 1, k + 2] < mean_max_df[iscen + 1, k + 2]) { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval2 <- t.test(fv2_frame_neutral[, k], alternative = "less", mu = curr_frame_a2[, k - 15])$p.value } else { # performing t tests of metric values from simulations with selection against # metric values from simulations without selection (neutral expectation) curr_pval2 <- t.test(fv2_frame_neutral[, k], alternative = "greater", mu = curr_frame_a2[, k - 15])$p.value } if (length(pval_list) >= k) { pval_list[[k]] <- c(pval_list[[k]], curr_pval) } else { pval_list[[k]] <- curr_pval } if (length(pval_list2) >= k) { pval_list2[[k]] <- c(pval_list2[[k]], curr_pval2) } else { pval_list2[[k]] <- curr_pval2 } if (curr_pval < 0.05) { sig_p_frame[ib, k] <- sig_p_frame[ib, k] + 1 } if (curr_pval2 < 0.05) { sig_p_frame[ib + 96, k] <- sig_p_frame[ib + 96, k] + 1 } } } # 8: FDRs FDR_list <- lapply(pval_list, function(x) p.adjust(x, method = "fdr")) FDR_list2 <- lapply(pval_list2, function(x) p.adjust(x, method = "fdr")) for (k in 16:18) { for (l in 1:length(FDR_list[[k]])) { if (!is.na(FDR_list[[k]][l])) { if (FDR_list[[k]][l] < 0.1) { sig_FDR_frame[ib, k] <- sig_FDR_frame[ib, k] + 1 } } } } for (k in 16:18) { for (l in 1:length(FDR_list2[[k]])) { if (!is.na(FDR_list2[[k]][l])) { if (FDR_list2[[k]][l] < 0.1) { sig_FDR_frame[ib + 96, k] <- sig_FDR_frame[ib + 96, k] + 1 } } } } glob_lengths[i] <- length(glob_vec_neutral) maxp1_freq_df_neutral[i, ] <- maxp1_df_neutral[i, ] maxp2_freq_df_neutral[i, ] <- maxp2_df_neutral[i, ] minp1_freq_df_neutral[i, ] <- minp1_df_neutral[i, ] minp2_freq_df_neutral[i, ] <- minp2_df_neutral[i, ] maxp1_df_neutral[i, ] <- maxp1_df_neutral[i, ] / 100 maxp2_df_neutral[i, ] <- maxp2_df_neutral[i, ] / 100 minp1_df_neutral[i, ] <- minp1_df_neutral[i, ] / 100 minp2_df_neutral[i, ] <- minp2_df_neutral[i, ] / 100 } # swapping the delta pi and H-scan columns (for consistency with figures in paper) sig_p_frame <- sig_p_frame[, c("V1", "V2", "V3", "V4", "V5", "V6", "V7", "V8", "V9", "V10", "V11", "V12", "V15", "V13", "V14", "V16", "V17", "V18")] sig_FDR_frame <- sig_FDR_frame[, c("V1", "V2", "V3", "V4", "V5", "V6", "V7", "V8", "V9", "V10", "V11", "V12", "V15", "V13", "V14", "V16", "V17", "V18")] ### sig_p_frame <- cbind(all_acc_df$scenario, sig_p_frame) sig_p_frame <- cbind(sig_p_frame, all_acc_df$recombination_sweep) sig_FDR_frame <- cbind(all_acc_df$scenario, sig_FDR_frame) sig_FDR_frame <- cbind(sig_FDR_frame, all_acc_df$recombination_sweep) colnames(sig_p_frame) <- c("scenario", colnames(fv_neutral)[1:18], "recombination_sweep") colnames(sig_FDR_frame) <- c("scenario", colnames(fv_neutral)[1:18], "recombination_sweep") # creating data frame with total number of simulations for each scenario total_numbers <- data.frame(matrix(rep(100, nrow(sig_p_frame) * ncol(sig_p_frame)), nrow = nrow(sig_p_frame), ncol = ncol(sig_p_frame))) # subtracting the number of NAs from the total number of simulations NA_frame <- cbind(sig_p_frame$scenario, NA_frame[, 1:18], sig_p_frame$recombination_sweep) total_numbers[, 2:19] <- total_numbers[, 2:19] - NA_frame # calculating proportions from the frequency values prop_p_frame <- sig_p_frame[, 2:19] / total_numbers[, 2:19] prop_FDR_frame <- sig_FDR_frame[, 2:19] / total_numbers[, 2:19] prop_p_frame <- cbind(sig_p_frame$scenario, prop_p_frame, sig_p_frame$recombination_sweep) prop_FDR_frame <- cbind(sig_p_frame$scenario, prop_FDR_frame, sig_p_frame$recombination_sweep) colnames(prop_p_frame) <- c("scenario", colnames(prop_p_frame)[2:19], "recombination_sweep") colnames(prop_FDR_frame) <- c("scenario", colnames(prop_FDR_frame)[2:19], "recombination_sweep") write.table(prop_p_frame, "rev_sig_sel2neutral_pval_18_proportions_central_v2.tsv", quote = FALSE, sep = "\t", row.names = FALSE) write.table(prop_FDR_frame, "rev_sig_sel2neutral_FDR_18_proportions_central_v2.tsv", quote = FALSE, sep = "\t", row.names = FALSE) write.table(sig_p_frame, "rev_sig_sel2neutral_pval_18_frequencies_central_v2.tsv", quote = FALSE, sep = "\t", row.names = FALSE) write.table(sig_FDR_frame, "rev_sig_sel2neutral_FDR_18_frequencies_central_v2.tsv", quote = FALSE, sep = "\t", row.names = FALSE) colnames(maxp1_freq_df_neutral) <- colnames(fv_neutral)[1:18] colnames(maxp2_freq_df_neutral) <- colnames(fv_neutral)[1:18] colnames(minp1_freq_df_neutral) <- colnames(fv_neutral)[1:18] colnames(minp2_freq_df_neutral) <- colnames(fv_neutral)[1:18] maxp1_freq_df_new_neutral <- data.frame("scenario" = title_vec_neutral, maxp1_freq_df_neutral) maxp2_freq_df_new_neutral <- data.frame("scenario" = title_vec_neutral, maxp2_freq_df_neutral) minp1_freq_df_new_neutral <- data.frame("scenario" = title_vec_neutral, minp1_freq_df_neutral) minp2_freq_df_new_neutral <- data.frame("scenario" = title_vec_neutral, minp2_freq_df_neutral) colnames(maxp1_df_neutral) <- colnames(fv_neutral)[1:18] colnames(maxp2_df_neutral) <- colnames(fv_neutral)[1:18] colnames(minp1_df_neutral) <- colnames(fv_neutral)[1:18] colnames(minp2_df_neutral) <- colnames(fv_neutral)[1:18] maxp1_df_new_neutral <- data.frame("scenario" = title_vec_neutral, maxp1_df_neutral) maxp2_df_new_neutral <- data.frame("scenario" = title_vec_neutral, maxp2_df_neutral) minp1_df_new_neutral <- data.frame("scenario" = title_vec_neutral, minp1_df_neutral) minp2_df_new_neutral <- data.frame("scenario" = title_vec_neutral, minp2_df_neutral) write.table(maxp1_df_new_neutral, "rev_18_prediction_freq_maxp1_1in1_neutral_v2.tsv", quote = FALSE, sep = "\t", row.names = FALSE) write.table(maxp2_df_new_neutral, "rev_18_prediction_freq_maxp2_1in1_neutral_v2.tsv", quote = FALSE, sep = "\t", row.names = FALSE) write.table(minp1_df_new_neutral, "rev_18_prediction_freq_minp1_1in1_neutral_v2.tsv", quote = FALSE, sep = "\t", row.names = FALSE) write.table(minp2_df_new_neutral, "rev_18_prediction_freq_minp2_1in1_neutral_v2.tsv", quote = FALSE, sep = "\t", row.names = FALSE) maxp1_freq_df_new_neutral <- data.frame(maxp1_freq_df_new_neutral, "recombination" = "high") maxp2_freq_df_new_neutral <- data.frame(maxp2_freq_df_new_neutral, "recombination" = "low") minp1_freq_df_new_neutral <- data.frame(minp1_freq_df_new_neutral, "recombination" = "high") minp2_freq_df_new_neutral <- data.frame(minp2_freq_df_new_neutral, "recombination" = "low") max_freq_df_neutral <- rbind(maxp1_freq_df_new_neutral, maxp2_freq_df_new_neutral) min_freq_df_neutral <- rbind(minp1_freq_df_new_neutral, minp2_freq_df_new_neutral) rownames(max_freq_df_neutral) <- 1:nrow(max_freq_df_neutral) rownames(min_freq_df_neutral) <- 1:nrow(min_freq_df_neutral) maxp1_df_new_neutral <- data.frame(maxp1_df_new_neutral, "recombination" = "high") maxp2_df_new_neutral <- data.frame(maxp2_df_new_neutral, "recombination" = "low") minp1_df_new_neutral <- data.frame(minp1_df_new_neutral, "recombination" = "high") minp2_df_new_neutral <- data.frame(minp2_df_new_neutral, "recombination" = "low") max_df_neutral <- rbind(maxp1_df_new_neutral, maxp2_df_new_neutral) min_df_neutral <- rbind(minp1_df_new_neutral, minp2_df_new_neutral) rownames(max_df_neutral) <- 1:nrow(max_df_neutral) rownames(min_df_neutral) <- 1:nrow(min_df_neutral) mean_max_df_neutral <- aggregate(max_df_neutral[, 2:19], list(max_df_neutral$recombination), mean) mean_min_df_neutral <- aggregate(min_df_neutral[, 2:19], list(min_df_neutral$recombination), mean) sd_max_df_neutral <- aggregate(max_df_neutral[, 2:19], list(max_df_neutral$recombination), sd) sd_min_df_neutral <- aggregate(min_df_neutral[, 2:19], list(min_df_neutral$recombination), sd) mean_acc_df_neutral <- mean_max_df_neutral sd_acc_df_neutral <- sd_max_df_neutral # determining if lower or higher values have higher predictive accuracy for (i in 1:nrow(mean_max_df)) { for (j in 3:ncol(mean_max_df)) { if (mean_min_df[i, j] > mean_max_df[i, j]) { mean_acc_df_neutral[i, j] <- mean_min_df_neutral[i, j] sd_acc_df_neutral[i, j] <- sd_min_df_neutral[i, j] } } } all_acc_df_neutral <- max_df_neutral all_freq_df_neutral <- max_freq_df_neutral for (i in 1:nrow(max_df)) { for (j in 2:(ncol(max_df) - 2)) { if (min_df[i, j] > max_df[i, j]) { all_acc_df_neutral[i, j] <- min_df_neutral[i, j] all_freq_df_neutral[i, j] <- min_freq_df_neutral[i, j] } } } colnames(mean_acc_df_neutral)[1] <- "recombination" colnames(sd_acc_df_neutral)[1] <- "recombination" mean_acc_df_neutral <- mean_acc_df_neutral[-c(3, 4), ] sd_acc_df_neutral <- sd_acc_df_neutral[-c(3, 4), ] write.table(mean_acc_df_neutral, "rev_mean_accuracy_18_prediction_1in1_neutral_v2.tsv", quote = FALSE, sep = "\t", row.names = FALSE) write.table(sd_acc_df_neutral, "rev_sd_accuracy_18_prediction_1in1_neutral_v2.tsv", quote = FALSE, sep = "\t", row.names = FALSE) all_acc_df_neutral <- all_acc_df_neutral[-c(97:192), ] write.table(all_acc_df_neutral, "rev_accuracy_18_prediction_1in1_neutral_v2.tsv", quote = FALSE, sep = "\t", row.names = FALSE) all_freq_df_neutral <- all_freq_df_neutral[-c(97:192), ] pval_df_neutral <- all_freq_df_neutral for (i in 1:nrow(all_freq_df_neutral)) { for (j in 2:(ncol(all_freq_df_neutral) - 1)) { if (i > 48) { i2 <- i - 48 } else { i2 <- i } curr_test <- prop.test(all_freq_df_neutral[i, j], 100, p = 0.05, alternative = "two.sided") pval_df_neutral[i, j] <- curr_test$p.value } } FDR_df_neutral <- pval_df_neutral pval_vec <- as.vector(unlist(pval_df_neutral[, -c(1, 30)])) FDR_vec <- p.adjust(pval_vec, "fdr") # all FDR values are larger than 0.1!!! no significant deviations from expectation!!! write.table(all_freq_df_neutral, "rev_18_prediction_1in1_freq_neutral_v2.tsv", quote = FALSE, sep = "\t", row.names = FALSE) pval_frame <- all_freq_df_neutral FDR_frame <- all_freq_df_neutral for (i in 1:nrow(pval_frame)) { for (j in 2:(ncol(pval_frame) - 2)) { pval_frame[i, j] <- pval_vec[((j - 1) * nrow(pval_frame)) + i] FDR_frame[i, j] <- FDR_vec[((j - 1) * nrow(pval_frame)) + i] } } write.table(pval_frame, "rev_18_prediction_1in1_freq_pval_neutral_v2.tsv", quote = FALSE, sep = "\t", row.names = FALSE) write.table(FDR_frame, "rev_18_prediction_1in1_freq_FDR_neutral_v2.tsv", quote = FALSE, sep = "\t", row.names = FALSE) ### ### end of neutral simulations part ### ### ### ### plots of accuracies derived from proportions of simulations that correctly identified sweep window ### ### ### heatmap for the above accuracies (1-in-1 window) all_acc_df <- data.frame(all_acc_df, "recombination" = str_detect(all_acc_df$recombination_sweep, "high", negate = FALSE)) all_acc_df$recombination <- str_replace(all_acc_df$recombination, "TRUE", "high") all_acc_df$recombination <- str_replace(all_acc_df$recombination, "FALSE", "low") all_acc_df <- data.frame(all_acc_df, "sweep" = str_detect(all_acc_df$recombination_sweep, "hard", negate = FALSE)) all_acc_df$sweep <- str_replace(all_acc_df$sweep, "TRUE", "hard") all_acc_df$sweep <- str_replace(all_acc_df$sweep, "FALSE", "soft") # reorder data frame (first by migration rate, second by recombination rate, third by growth ("population_size"), fourth by demography, fifth by sweep type, sixth by "equal_size") new_all_acc_df <- all_acc_df[-which(all_acc_df$recombination == "low"), ] new_all_acc_df <- new_all_acc_df[-which(new_all_acc_df$equal_size == "unequal"), ] new_all_acc_df <- new_all_acc_df[with(new_all_acc_df, order(migration_rate, sweep, population_size, demography)), ] pal_breaks <- seq(0, 1, 1/100) png("rev_accuracy_heatmap_scale_byr_16_prediction_1in1_wo_lowrec_unequal_order_mig_sweep_growth_dem_v2.png", width = 1300, height = 1300) all_acc_mat <- as.matrix(new_all_acc_df[, -c(1, 18:ncol(new_all_acc_df))]) rownames(all_acc_mat) <- new_all_acc_df$migration_rate colnames(all_acc_mat) <- colnames(new_all_acc_df)[-c(1, 18:ncol(new_all_acc_df))] pheatmap::pheatmap(all_acc_mat, scale = "none", cluster_rows = FALSE, cluster_cols = FALSE, treeheight_row = 0, treeheight_col = 0, breaks = pal_breaks, fontsize = 15, color = colorRampPalette(rev(brewer.pal(n = 11, name ="RdYlBu")))(100)) dev.off() ### new_all_acc_df <- all_acc_df[-which(all_acc_df$recombination == "high"), ] new_all_acc_df <- new_all_acc_df[-which(new_all_acc_df$equal_size == "unequal"), ] new_all_acc_df <- new_all_acc_df[with(new_all_acc_df, order(migration_rate, sweep, population_size, demography)), ] png("rev_accuracy_heatmap_scale_byr_16_prediction_1in1_wo_highrec_unequal_order_mig_sweep_growth_dem_v2.png", width = 1300, height = 1300) all_acc_mat <- as.matrix(new_all_acc_df[, -c(1, 18:ncol(new_all_acc_df))]) rownames(all_acc_mat) <- new_all_acc_df$migration_rate colnames(all_acc_mat) <- colnames(new_all_acc_df)[-c(1, 18:ncol(new_all_acc_df))] pheatmap::pheatmap(all_acc_mat, scale = "none", cluster_rows = FALSE, cluster_cols = FALSE, treeheight_row = 0, treeheight_col = 0, breaks = pal_breaks, fontsize = 15, color = colorRampPalette(rev(brewer.pal(n = 11, name ="RdYlBu")))(100)) dev.off() ### new_all_acc_df <- all_acc_df[-which(all_acc_df$recombination == "low"), ] new_all_acc_df <- new_all_acc_df[-which(new_all_acc_df$equal_size == "equal"), ] new_all_acc_df <- new_all_acc_df[with(new_all_acc_df, order(migration_rate, sweep, population_size, demography)), ] png("rev_accuracy_heatmap_scale_byr_16_prediction_1in1_wo_lowrec_equal_order_mig_sweep_growth_dem_v2.png", width = 1300, height = 1300) all_acc_mat <- as.matrix(new_all_acc_df[, -c(1, 18:ncol(new_all_acc_df))]) rownames(all_acc_mat) <- new_all_acc_df$migration_rate colnames(all_acc_mat) <- colnames(new_all_acc_df)[-c(1, 18:ncol(new_all_acc_df))] pheatmap::pheatmap(all_acc_mat, scale = "none", cluster_rows = FALSE, cluster_cols = FALSE, treeheight_row = 0, treeheight_col = 0, breaks = pal_breaks, fontsize = 15, color = colorRampPalette(rev(brewer.pal(n = 11, name ="RdYlBu")))(100)) dev.off() ### new_all_acc_df <- all_acc_df[-which(all_acc_df$recombination == "high"), ] new_all_acc_df <- new_all_acc_df[-which(new_all_acc_df$equal_size == "equal"), ] new_all_acc_df <- new_all_acc_df[with(new_all_acc_df, order(migration_rate, sweep, population_size, demography)), ] png("rev_accuracy_heatmap_scale_byr_16_prediction_1in1_wo_highrec_equal_order_mig_sweep_growth_dem_v2.png", width = 1300, height = 1300) all_acc_mat <- as.matrix(new_all_acc_df[, -c(1, 18:ncol(new_all_acc_df))]) rownames(all_acc_mat) <- new_all_acc_df$migration_rate colnames(all_acc_mat) <- colnames(new_all_acc_df)[-c(1, 18:ncol(new_all_acc_df))] pheatmap::pheatmap(all_acc_mat, scale = "none", cluster_rows = FALSE, cluster_cols = FALSE, treeheight_row = 0, treeheight_col = 0, breaks = pal_breaks, fontsize = 15, color = colorRampPalette(rev(brewer.pal(n = 11, name ="RdYlBu")))(100)) dev.off() ### ### accuracy calculations for 1 in 3 windows ### # parsing input file names filename_vec <- sapply(title_vec, function(x) gsub(" ", "_", x)) filename_vec <- sapply(filename_vec, function(x) str_remove(x, "\\(")) filename_vec <- sapply(filename_vec, function(x) str_remove(x, "\\)")) filename_vec <- sapply(filename_vec, function(x) str_remove(x, ",")) filename_vec <- sapply(filename_vec, function(x) str_remove(x, ";")) filename_vec <- sapply(filename_vec, function(x) gsub("0\\.", "0_", x)) filename_vec_neutral <- sapply(title_vec_neutral, function(x) gsub(" ", "_", x)) filename_vec_neutral <- sapply(filename_vec_neutral, function(x) str_remove(x, "\\(")) filename_vec_neutral <- sapply(filename_vec_neutral, function(x) str_remove(x, "\\)")) filename_vec_neutral <- sapply(filename_vec_neutral, function(x) str_remove(x, ",")) filename_vec_neutral <- sapply(filename_vec_neutral, function(x) str_remove(x, ";")) filename_vec_neutral <- sapply(filename_vec_neutral, function(x) gsub("0\\.", "0_", x)) # initialising data frames maxp1_df <- as.data.frame(matrix(rep(0, 18 * length(title_vec)), length(title_vec), 18)) maxp2_df <- as.data.frame(matrix(rep(0, 18 * length(title_vec)), length(title_vec), 18)) minp1_df <- as.data.frame(matrix(rep(0, 18 * length(title_vec)), length(title_vec), 18)) minp2_df <- as.data.frame(matrix(rep(0, 18 * length(title_vec)), length(title_vec), 18)) # looping through all scenarios for (i in 1:length(title_vec)) { print(i) glob_vec <- Sys.glob(file.path(paste(scenario_vec[i]))) glob_vec_p3 <- Sys.glob(file.path(paste(scenario_vec_p3[i]))) glob_vec2 <- Sys.glob(file.path(paste(path, fst_vec[i], sep = ""))) glob_vec3 <- Sys.glob(file.path(paste(path2, beta_vec[i], sep = ""))) glob_vec4 <- Sys.glob(file.path(paste(path, hscan_vec[i], ".windows", sep = ""))) df_list <- list() for (j in 1:100) { # reading in metric input files fv <- read.table(glob_vec[j], header = TRUE) fv_p3 <- read.table(glob_vec_p3[j], header = TRUE) statvars <- colsplit(colnames(fv), "_win", c("metric", "window")) fv <- unlist(sapply(fv, unlist)) fv_p3 <- unlist(sapply(fv_p3, unlist)) fv <- as.data.frame(matrix(fv, nrow = 50)) fv_p3 <- as.data.frame(matrix(fv_p3, nrow = 50)) rownames(fv) <- 1:50 rownames(fv_p3) <- 1:50 colnames(fv) <- unique(statvars[, 1]) colnames(fv_p3) <- unique(statvars[, 1]) hscan_frame <- read.table(glob_vec4[j], header = TRUE) fv <- data.frame(fv, "hscan" = hscan_frame[1:50, 2]) fv <- data.frame(fv, "delta_pi" = fv[1:50, 1] - fv_p3[1:50, 1]) fst_frame <- read.table(glob_vec2[j], header = TRUE) fv <- data.frame(fv, "Fst" = fst_frame[1:50, 6]) beta_frame <- readRDS(glob_vec3[j]) rownames(beta_frame) <- 1:50 fv <- data.frame(fv, "beta_q0" = beta_frame[1:50, 2]) fv <- data.frame(fv, "beta_q1" = beta_frame[1:50, 3]) fv <- data.frame(fv, "beta_q2" = beta_frame[1:50, 4]) fv <- data.frame(fv, "window_num" = (1:50) * 4000050 / 50) fv1 <- fv[1:25, ] fv2 <- fv[26:50, ] for (k in 1:(ncol(fv) - 1)) { part1 <- fv1[, k] part2 <- fv2[, k] maxp1_1 <- order(part1, decreasing = TRUE)[1] maxp1_2 <- order(part1, decreasing = TRUE)[2] maxp2_1 <- order(part2, decreasing = TRUE)[1] maxp2_2 <- order(part2, decreasing = TRUE)[2] minp1_1 <- order(part1, decreasing = FALSE)[1] minp1_2 <- order(part1, decreasing = FALSE)[2] minp2_1 <- order(part2, decreasing = FALSE)[1] minp2_2 <- order(part2, decreasing = FALSE)[2] # determining maxima and minima of metrics if (maxp1_1 %in% 12:14) { maxp1_df[i, k] <- maxp1_df[i, k] + 1 } if (maxp2_1 %in% 12:14) { maxp2_df[i, k] <- maxp2_df[i, k] + 1 } if (minp1_1 %in% 12:14) { minp1_df[i, k] <- minp1_df[i, k] + 1 } if (minp2_1 %in% 12:14) { minp2_df[i, k] <- minp2_df[i, k] + 1 } } } maxp1_df[i, ] <- maxp1_df[i, ] / length(glob_vec) maxp2_df[i, ] <- maxp2_df[i, ] / length(glob_vec) minp1_df[i, ] <- minp1_df[i, ] / length(glob_vec) minp2_df[i, ] <- minp2_df[i, ] / length(glob_vec) } colnames(maxp1_df) <- colnames(fv)[1:18] colnames(maxp2_df) <- colnames(fv)[1:18] colnames(minp1_df) <- colnames(fv)[1:18] colnames(minp2_df) <- colnames(fv)[1:18] maxp1_df_new <- data.frame("scenario" = title_vec, maxp1_df) maxp2_df_new <- data.frame("scenario" = title_vec, maxp2_df) minp1_df_new <- data.frame("scenario" = title_vec, minp1_df) minp2_df_new <- data.frame("scenario" = title_vec, minp2_df) write.table(maxp1_df_new, "rev_18_prediction_freq_maxp1_1in3.tsv", quote = FALSE, sep = "\t", row.names = FALSE) write.table(maxp2_df_new, "rev_18_prediction_freq_maxp2_1in3.tsv", quote = FALSE, sep = "\t", row.names = FALSE) write.table(minp1_df_new, "rev_18_prediction_freq_minp1_1in3.tsv", quote = FALSE, sep = "\t", row.names = FALSE) write.table(minp2_df_new, "rev_18_prediction_freq_minp2_1in3.tsv", quote = FALSE, sep = "\t", row.names = FALSE) maxp1_df_new <- data.frame(maxp1_df_new, "recombination" = "high") maxp2_df_new <- data.frame(maxp2_df_new, "recombination" = "low") minp1_df_new <- data.frame(minp1_df_new, "recombination" = "high") minp2_df_new <- data.frame(minp2_df_new, "recombination" = "low") max_df <- rbind(maxp1_df_new, maxp2_df_new) min_df <- rbind(minp1_df_new, minp2_df_new) rownames(max_df) <- 1:nrow(max_df) rownames(min_df) <- 1:nrow(min_df) max_df <- data.frame(max_df, "hard_sweep" = str_detect(max_df$scenario, "hard", negate = FALSE)) min_df <- data.frame(min_df, "hard_sweep" = str_detect(min_df$scenario, "hard", negate = FALSE)) mean_max_df <- aggregate(max_df[, 2:19], list(max_df$recombination, max_df$hard_sweep), mean) mean_min_df <- aggregate(min_df[, 2:19], list(min_df$recombination, min_df$hard_sweep), mean) sd_max_df <- aggregate(max_df[, 2:19], list(max_df$recombination, max_df$hard_sweep), sd) sd_min_df <- aggregate(min_df[, 2:19], list(min_df$recombination, min_df$hard_sweep), sd) colnames(mean_max_df)[1:2] <- c("recombination", "hard_sweep") colnames(mean_min_df)[1:2] <- c("recombination", "hard_sweep") colnames(sd_max_df)[1:2] <- c("recombination", "hard_sweep") colnames(sd_min_df)[1:2] <- c("recombination", "hard_sweep") mean_acc_df <- mean_max_df sd_acc_df <- sd_max_df mean_acc_df <- data.frame(mean_acc_df, "recombination_sweep" = paste(mean_acc_df$recombination, mean_acc_df$hard_sweep, sep = "_")) sd_acc_df <- data.frame(sd_acc_df, "recombination_sweep" = paste(sd_acc_df$recombination, sd_acc_df$hard_sweep, sep = "_")) mean_acc_df$recombination_sweep <- str_replace(mean_acc_df$recombination_sweep, "FALSE", "soft") mean_acc_df$recombination_sweep <- str_replace(mean_acc_df$recombination_sweep, "TRUE", "hard") sd_acc_df$recombination_sweep <- str_replace(sd_acc_df$recombination_sweep, "FALSE", "soft") sd_acc_df$recombination_sweep <- str_replace(sd_acc_df$recombination_sweep, "TRUE", "hard") # determining if lower or higher values have higher predictive accuracy for (i in 1:nrow(mean_max_df)) { for (j in 3:ncol(mean_max_df)) { if (mean_min_df[i, j] > mean_max_df[i, j]) { mean_acc_df[i, j] <- mean_min_df[i, j] sd_acc_df[i, j] <- sd_min_df[i, j] } } } all_acc_df <- max_df all_acc_df <- data.frame(all_acc_df, "recombination_sweep" = paste(all_acc_df$recombination, all_acc_df$hard_sweep, sep = "_")) all_acc_df$recombination_sweep <- str_replace(all_acc_df$recombination_sweep, "FALSE", "soft") all_acc_df$recombination_sweep <- str_replace(all_acc_df$recombination_sweep, "TRUE", "hard") for (i in 1:nrow(max_df)) { for (j in 2:(ncol(max_df) - 2)) { if (min_df[i, j] > max_df[i, j]) { all_acc_df[i, j] <- min_df[i, j] } } } mean_acc_df <- mean_acc_df[, -c(1, 2)] sd_acc_df <- sd_acc_df[, -c(1, 2)] write.table(mean_acc_df, "rev_mean_accuracy_18_prediction_1in3.tsv", quote = FALSE, sep = "\t", row.names = FALSE) write.table(sd_acc_df, "rev_sd_accuracy_18_prediction_1in3.tsv", quote = FALSE, sep = "\t", row.names = FALSE) write.table(all_acc_df, "rev_accuracy_18_prediction_1in3.tsv", quote = FALSE, sep = "\t", row.names = FALSE) ### ###