In eriqande/inbredPermute: Simple permutation procedures for distributions of rxy values amongst parents

inbredPermute

Simple scripts/procedures in development to permute values of rxy around to see if survival is related to inbreeding in salmon pedigrees. Joint work with David Vendrami.

Basic workflow

The basic idea here is that we should be able to proceed directly with two ingredients:

1. A raw kingroup output file that includes $r_{xy}$ values. The top of that file should look something like this: r cat(readLines(system.file("data_files/MKHM_kingroup_output.csv", package="inbredPermute", mustWork = T), n = 8), sep="\n")

2. A list of trios that were found in the population, that include offspring whose parents are the ones that were included in the kingroup analysis. I must have three columns that are named Ma, Pa, and Kid (which are case-sensitive), and it can have any other columns that you might happen to have in there, too. Here is an example: r cat(readLines(system.file("data_files/MKHM_trios.txt", package="inbredPermute", mustWork = T), n = 10), sep="\n")

3. A file with spawn date and sex of individuals in this sort of format: r cat(readLines(system.file("data_files/MKHM2011_metadata.txt", package="inbredPermute", mustWork = T), n = 10), sep="\n")

Then, to run the function you would do something like:

    library(inbredPermute)
    do_perm_output <- do_the_perms(kg_file = system.file("data_files/MKHM_kingroup_output.csv", package="inbredPermute", mustWork = T),
                           trio_file = system.file("data_files/MKHM_trios.txt", package="inbredPermute", mustWork = T),
                           meta_file = system.file("data_files/MKHM2011_metadata.txt", package="inbredPermute", mustWork = T),
                           ItalianCommas = TRUE,
                           REPS = 1000  # you will want to do more in practice
                           )

    do_perm_output[1:3]
    head(do_perm_output$Obs)
    do_perm_output$Simmed[1:10, 1:10]

and here were are working on comparing each simulated distribution to the observed distribution of r_xy values:

# pass in a simulated vector of r_xy called y.
# This does a Mann_Whitney U-test on those dsns.
wilx_func <- function(y) unlist(wilcox.test(as.numeric(do_perm_output$Obs$mapa_rxy), y, alternative = "less")[c("p.value", "statistic")])

wilx_tests <- t(apply(do_perm_output$Simmed, 2, wilx_func))
hist(wilx_tests[, 1])  # clearly not much going on there.

# what if, instead we asked about the max inbreeding coeff?
max_obs <- max(as.numeric(do_perm_output$Obs$mapa_rxy))
# number of simmed inds with higher inbreeding than any observed
gt_obs <- apply(do_perm_output$Simmed, 2, function(x) sum(x > max_obs))
hist(gt_obs)
mean(gt_obs == 0)  # this is a sort of "p-value"  

# finally, what if we look at the number of offspring with inbreeding coeff 
# greater than 0.25?
num_obs_gt_.25 <- sum(as.numeric(do_perm_output$Obs$mapa_rxy) > 0.25)
sim_obs_gt_.25 <- apply(do_perm_output$Simmed, 2, function(x) sum(x > 0.25))

Terms

As a work partially of the United States Government, this package is in the public domain within the United States.

Of course, it would be awfully lame of anyone to take anything found within here and use it as their own before we tried publishing any of this, should we choose to do that.

See TERMS.md for more information.

eriqande/inbredPermute documentation built on May 16, 2019, 8:47 a.m.