Utilities for Making NewHybs input Files

Intro/Overview

After I write some functions for segregating haplotypes and creating hybrids I am going to need some functions for turning data frames of genotypes into newhybrids input files. I am going to work on this before doing the segregation code since I think it will inform what I want the end result of that segregation code to look like.

Our basic goal here is to pass in a data frame of individuals that has at least the columns: id, variant, gene_copy, allele, and also pass in a vector of variant names, and then this will find all the common variant names and make a data frame that is suitable for printing out for newhybrids.

Mostly this will be a lot of tidyr…

library(tidyverse)

Loading tidyverse: ggplot2
Loading tidyverse: tibble
Loading tidyverse: tidyr
Loading tidyverse: readr
Loading tidyverse: purrr
Loading tidyverse: dplyr
Conflicts with tidy packages --------------------------------------------------------------------------------
filter(): dplyr, stats
lag():    dplyr, stats

library(stringr)
source("R/nh_tablify.R")
source("R/write_nh.R")

Functions

This is just one that makes tables of nonLumped genotypes that we will be able to print to newhybrids files. (NOTE to self, this is going to work poorly if the group is missing data at all indivs at a locus you want…)

#' @param D a tidy data frame with at least the columns: id, variant, gene_copy, allele
#' @param V a vector of variant names giving which variants to use and the order they are
#' to be in.
#' @param opt_str  The options string to be added to each individual (like "z0s")
#' @param id_prepend this string will be prepended to the fish id (i.e. use it for "train_", or "test_")
nh_tablify <- function(D, V, opt_str = "", id_prepend = "") {
  D1 <- D %>%
    filter(variant %in% V)
  
  common_loc <- V[V %in% D1$variant]  # these are the loci from V that are in D, in V-order
  
  # now spread the gene copies and then the variants
  D1 %>%
    select(id, variant, gene_copy, allele) %>%
    tidyr::spread(data = ., key = gene_copy, value = allele) %>%
    mutate(gc1 = as.integer(ifelse(is.na(`1`), -1, `1`)),
           gc2 = as.integer(ifelse(is.na(`2`), -1, `2`))) %>%
    select(-`1`, -`2`) %>%
    mutate(geno = paste(gc1, gc2, sep = " ")) %>%
    select(-gc1, -gc2) %>% 
    tidyr::spread(data = ., key = variant, value = geno, fill = "-1 -1") %>%
    .[, c("id", common_loc)] %>%
    mutate(option = opt_str) %>%
    mutate(id = paste(id_prepend, id, sep = "")) %>%
    select(id, option, everything())
  
}

Now, if you have bind_rows-ed together a bunch of products of nh_tablify() then you can pass them to this function to write a newhybs file:

#' write a new_hybs file
#'
#' @param G a data frame like that prepared by nh_tablify, or perhaps several of those
#' that have been bind_rows()-ed together.
#' @param path  the file path you want to write to.
write_nh <- function(G, path = "nh_data.txt") {
  # first write the preamble
  cat("NumIndivs ", nrow(G), "\nNumLoci ", ncol(G) - 2, "\nDigits 1\nFormat NonLumped\n\n", file = path)
  
  # then the locus names
  cat(c("LocusNames  ", names(G)[-(1:2)]), sep = " ", eol = "\n", file = path, append = TRUE)
  
  # then format the genotype frame and write it
  G %>%
    mutate(id = paste("n", id, sep = " ")) %>%
    ungroup() %>%
    mutate(idx = 1:n()) %>%
    select(idx, everything()) %>%
    write.table(file = path, append = TRUE, quote = FALSE, row.names = FALSE, col.names = FALSE, sep = "   ", eol = "\n")
  
}

Putting those together

OK, now that we have those functions we can put them all together with a simple example using the intermediate data stored from the last document (02). First let us get the data.

TT <- readRDS("intermediates/02/test_and_train.rds")
MM <- readRDS("intermediates/02/marker_rankings.rds")

Now, choose the top 1000 loci to use:

v1000 <- MM %>%
  filter(selectable == TRUE & cumsum(selectable) <= 1000) %>%
  .$variant

Now, we are going are just going to put pure wild and pure farmed individuals in there as a test case. We are going to include the training individuals with the z0s and z1s options. We can do this all into one big list first:

nh_list <- list(
train_farm = TT %>%
  filter(test_or_train == "train", group == "farmed") %>%
  nh_tablify(., v1000, opt_str = "z0s", id_prepend = "train_"),
train_wild = TT %>%
  filter(test_or_train == "train", group == "wild") %>%
  nh_tablify(., v1000, opt_str = "z1s", id_prepend = "train_"),
test_farm = TT %>%
  filter(test_or_train == "test", group == "farmed") %>%
  nh_tablify(., v1000, opt_str = "", id_prepend = "test_"),
test_wild = TT %>%
  filter(test_or_train == "test", group == "wild") %>%
  nh_tablify(., v1000, opt_str = "", id_prepend = "test_")
)
dir.create("scratch", showWarnings = FALSE)
bind_rows(nh_list) %>%
  write_nh(path = "scratch/all_pures.txt")

Note that if I want to run that in newhybrids with a fixed prior for the different categories, I can do that like this:

~/Documents/git-repos/newhybrids/newhybs -d all_pures.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

Looking at the output, it is easy to verify that when the prior is fixed and equal for all categories the aa-PofZ.txt and aa-ScaledLikelihood.txt files contain exactly the same values. But the aa-PofZ.txt also has an IndivName column (using the “n” option for names), so that will be the one that we want to use.

Great. We are almost ready to put it all together, I just need to make a function that creates hybrids from the test individuals.