Making Multiple Data Sets

Intro/Overview

Just going to show how to create lots of data sets that can then be run through new hybrids. We create our own nomenclature so that the directory name tells us about the simulation.

library(tidyverse)
library(stringr)
library(SalarHybPower)

We also want to have some example data. We will just use the Nova Scotia data. But we will call this dat:

dat <- readRDS("intermediates/01/tidy-west.rds") %>%
  mutate(pop = str_sub(id, 1, 3),
         group = ifelse(pop == "AQU" | pop == "WLN", "farmed", "wild"))

First, count how many individuals we have from each wild population:

dat %>%
  group_by(pop) %>%
  summarize(num = n_distinct(id))

This means that if we are doing F2’s for LHR we will only be able to make 4 or 5 per time, and for backcrosses, we will only get about 4 per data set out of them. OK. Well, let us not worry about the number of individuals we get from each population. Let’s just do SPLITS different splits of the data into training and holdout…

Set up some info for the sims:

SPLITS <- 3 # the number of times to split the data and rank markers
LOCS <- c(24, 48, 144, 192, 480, 720, 1000)  # number of loci
POPLIST <- c("BDN", "CNR", "GRR", "LHR", "LPR")  # names of wild populations
HYB_CATS <- c("PureW", "PureF", "F1", "F2", "BX")

Then just cycle over things and do it

main_out_dir <- "nh_reps_directory"
dir.create(main_out_dir)
set.seed(555)
for (s in 1:SPLITS) {
  SAR <- split_and_rank(dat)
  for (pop in POPLIST) {
    for (hc in HYB_CATS) {
      for (locs in LOCS) {
        dirname <- paste(s,pop,hc,locs, sep = "_")
        create_hybrid_dataset(SAR = SAR, wild_pop = pop, hyb_cat = hc, L = locs, 
                              dir = file.path(main_out_dir, dirname))
      }
    }
  }
  
}

That takes a few minutes and creates a large number of directories, (like 525 or so), each with a single simulated data set in it.
The directories are named like this: 1_GRR_PureF_480 which denotes:

• Split = 1
• Population = GRR
• Simulated hybrid category = Pure Farmed
• Number of loci = 480

Now I just need to run NewHybrids over each like this:

~/Documents/git-repos/newhybrids/newhybs -d nh_data.txt -g P0 1 0 0 -g p1 0 0 1 -g F1 0 1 0 -g F2 .25 .5 .25 -g BX0 .5 .5 0 -g BX1 0 .5 .5 --pi-prior fixed  1 1 1 1 1 1

We will do that using the GNU parallel script.

Creating a GNU parallel script

We could use the parallel R package, but I have had better success using GNU parallel – a Perl script. I have included it in the bin directory in the repo.

We basically need to write a series of shell commands, each one launching newhybrids on a different data set.

We can use R to write those out. Note that this all assumes that the nh_reps_directory is in the top level of the repository, and that the command to run parallel will be given in that nh_reps_directory. Note that it is imporant to divert stderr to a newhybs_stderr file so we can search for cases that had underflow issues. I am going to do 100 burn in and 500 sweeps, cuz we can do short runs when we are working with having some indivs of known origin.

comms <- lapply(dir("nh_reps_directory"), function(x) {
  paste0("echo \"Starting ", 
         x, 
         " at $(date)\"; cd ", 
         x, 
         "; ../../bin/newhybs -d nh_data.txt -g P0 1 0 0 -g p1 0 0 1 -g F1 0 1 0 -g F2 .25 .5 .25 -g BX0 .5 .5 0 -g BX1 0 .5 .5 --pi-prior fixed  1 1 1 1 1 1 --no-gui --burn-in 100 --num-sweeps 500 --seeds ",
         paste(ceiling(runif(2, min = 100, max = 10000000)), collapse = " "), 
         " > newhybs_stdout.txt 2> newhybs_stderr.txt; cd ..; echo \"Done with ", x, " at $(date)\"")
})
cat(unlist(comms), sep = "\n", file = "para-comms.txt")

I synced all that to our 24 core box and then put the runs on 22 cores:

2017-06-20 14:44 /nh_reps_directory/--% (master) pwd
/Users/eriq/Documents/git-repos/SalarHybPower/nh_reps_directory

# then put it on 22 cores:
2017-06-20 14:45 /nh_reps_directory/--% (master) cat ../para-comms.txt | ../bin/parallel -P22 > ../BIG_LOG.txt &

When doing such short runs, this seems to take less than an hour…

Slurping up the output

Since this is on the big machine at work, it will be easiest for me to grab the output using awk over ssh:

2017-06-21 05:57 /nh_reps_directory/--% (master) pwd
/Users/eriq/Documents/git-repos/SalarHybPower/nh_reps_directory

for i in *; do 
  warn=$(wc -l $i/newhybs_stderr.txt | awk '{print $1}'); 
  awk -v w=$warn -v i=$i 'NR>1 && !/train_/ {print i, w, $0}' $i/aa-PofZ.txt;
done  > ../3_splits_output.txt

Then I gzipped that and brought it over to the outputs folder in the repository (but which I won’t commit…)

Reading in the output, and making a quick plot

First, a job for readr:

nh_output <- read_table2("outputs/3_splits_output.txt.gz",   # requires readr 1.1.1
                         col_names = c("sim", "warn", "idx", "iname", 
                                       "pure_farmed", "pure_wild", "F1", "F2", "bx_farm", "bx_wild"))

Parsed with column specification:
cols(
  sim = col_character(),
  warn = col_integer(),
  idx = col_integer(),
  iname = col_character(),
  pure_farmed = col_double(),
  pure_wild = col_double(),
  F1 = col_double(),
  F2 = col_double(),
  bx_farm = col_double(),
  bx_wild = col_double()
)

Then separate some columns and then gather

nh_sepped <- nh_output %>%
  separate(col = sim, into = c("split", "pop", "true_hyb_cat", "num_loci"), sep = "_", convert = TRUE) %>%
  separate(col = iname, into = c("drop", "s_idx"), sep = "_", convert = TRUE) %>% 
  mutate(s_idx = paste("split", split, s_idx, sep = "_")) %>%
  select(-drop) 

Really quickly, check which runs received warnings:

nh_sepped %>%
  filter(warn > 0)

OK, all of them have num_loci = 1000.

Let’s just filter out the ones that had warnings and gather it, and make a infferred_hyb_cat a factor so it comes out in a good order in the plots.

nh_tidy <- nh_sepped %>% 
  filter(warn == 0) %>% 
  gather(key = "inferred_hyb_cat", value = "post_prob", pure_farmed:bx_wild) %>%
  mutate(inferred_hyb_cat = factor(inferred_hyb_cat, levels = c("bx_wild", "bx_farm", "F2", "F1", "pure_wild", "pure_farmed")))
# and save this
saveRDS(nh_tidy, file = "outputs/3_splits_nh_tidy.rds", compress = "xz")

And now we should be able to plot it.

# some colors that might work out better than ggplots defaults
hc_colors <- c(
  pure_farmed = rgb(152, 78, 163, maxColorValue = 255), # purple
  pure_wild = rgb(77, 175, 74, maxColorValue = 255), # green
  F1 = rgb(255, 255, 51, maxColorValue = 255),  # yellow
  F2 = rgb(247, 129, 191, maxColorValue = 255), # pink
  bx_wild = rgb(255, 127, 0, maxColorValue = 255),  # orange
  bx_farm = rgb(228, 26, 28, maxColorValue = 255) # red
)
nh_list <- nh_tidy %>% 
  split(., .$true_hyb_cat)
gg_list <- lapply(nh_list, function(x) {
  g <- ggplot(x, aes(x = s_idx, y = post_prob, fill = inferred_hyb_cat)) +
    geom_col() +
    facet_grid(num_loci ~ pop) +
    scale_fill_manual(values = hc_colors)
  
  g
})
dump <- lapply(names(gg_list), function(n) {
  ggsave(gg_list[[n]], 
         filename = paste0("outputs/first_look_", n, ".pdf"),
         width = 10, 
         height = 10)})

I will send those figures to Brendan.

From my first looking over of them, there are a couple of trends it seems to me:

• Populations with fewer individuals have fewer simulated individuals—that is why there are gaps in there. (though some of the big gaps at 1000 loci are due to the fact that newhybrids has underflow problems with lots of loci.)

• It is important to understand that I simulated individuals from separate populations rather than slurrying their alleles together to make a mean-wild-allele-freq pool. From the results, it is clear that different populations must have different allele frequencies, such that when all the training individuals are used as reference wild and farmed fish with the z option in NewHybrids, some wild fish of some populations end up looking like bx_wilds, etc. and some pure Farmed fish look like bx_farms. The latter is not so bad, I would think, because you don’t really expect bx_farmed to occur, as often in the populations.

• It is interesting that the truly-pure farmed fish are incorrectly inferred to be bx_farm fish at a higher rate with more loci. I have checked my simulation programs and I don’t think they are incorrect. I think it might have to do with using a mixture of populations in the fish with the z option.

• It might be better to run newhybrids with wild training individuals being included and getting the z option only if they are from the same population as the sampled (simulated) hybrids. This is the next thing I need to get on, but wanted to send these results on for now.

• For some of the populations, like BDN, CNR, and LHR it looks like hybrids/farmed fish can be identified pretty reliably. For GRR and LPR, it is a little tougher as there is a lot of overlap between pure wild and backcrossed fish.