context("Loading data") # Pedigree 1: # # G3 11 # / \ # G2 9 10 # / \ | # G1 6 7 8 # | /\ /\ # G0 1 2 3 4 5 # # # Pedigree 2: # # G3 39 _______ # / \ \ # G1 36 37 38 # |\ / | \ # G0 31 32 33 34 35 loaded_data <- structure(list(pid = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 31L, 32L, 33L, 34L, 35L, 36L, 37L, 38L, 39L), pid_dad = c(6L, 7L, 7L, 8L, 8L, 9L, 9L, 10L, 11L, 11L, 0L, 36L, 36L, 37L, 37L, 37L, 39L, 39L, 39L, 0L)), class = c("tbl_df", "tbl", "data.frame"), row.names = c(NA, -20L)) test_pop <- load_individuals(pid = loaded_data$pid, pid_dad = loaded_data$pid_dad, progress = FALSE) test_that("test_create_population works", { expect_failure(expect_null(test_pop)) expect_output(print(test_pop), regexp = "^Population with 20 individuals$") expect_equal(pop_size(test_pop), 20L) }) indvs <- get_individuals(test_pop) test_that("get_individuals works", { expect_failure(expect_null(indvs)) expect_equal(length(indvs), 20L) }) test_that("generation assignment works", { # 0, last generation: expect_equal(0L, unique(unlist(lapply(1L:5L, function(pid) get_generation(get_individual(test_pop, pid))))) ) expect_equal(0L, unique(unlist(lapply(31L:35L, function(pid) get_generation(get_individual(test_pop, pid))))) ) # 1, second last generation: expect_equal(1L, unique(unlist(lapply(6L:8L, function(pid) get_generation(get_individual(test_pop, pid))))) ) expect_equal(1L, unique(unlist(lapply(36L:38L, function(pid) get_generation(get_individual(test_pop, pid))))) ) # 2, third last generation: expect_equal(2L, unique(unlist(lapply(9L:10L, function(pid) get_generation(get_individual(test_pop, pid))))) ) expect_equal(2L, unique(unlist(lapply(39L, function(pid) get_generation(get_individual(test_pop, pid))))) ) # 3, fourth last generation: expect_equal(3L, unique(unlist(lapply(11L, function(pid) get_generation(get_individual(test_pop, pid))))) ) }) peds <- build_pedigrees(test_pop, progress = FALSE) test_that("build_pedigrees works", { expect_output(print(peds), regexp = "^List of 2 pedigrees \\(of size 11, 9\\)$") expect_equal(pedigrees_count(peds), 2L) }) ped <- peds[[1L]] pids <- sort(get_pids_in_pedigree(ped)) test_that("pedigree pids works", { expect_equal(length(pids), 11L) expect_true(all(pids == 1L:11L)) }) # ped 1: pids <- sort(get_pids_in_pedigree(peds[[1L]])) test_that("pedigree pids works", { expect_equal(length(pids), 11L) expect_true(all(pids == 1L:11L)) }) # ped 2: pids <- sort(get_pids_in_pedigree(peds[[2L]])) test_that("pedigree pids works", { expect_equal(length(pids), 9L) expect_true(all(pids == 31L:39L)) }) test_that("meiotic_dist works", { if (FALSE) { plot(peds) } expect_equal(4L, meiotic_dist(get_individual(test_pop, pid = 1L), get_individual(test_pop, pid = 10L))) expect_equal(6L, meiotic_dist(get_individual(test_pop, pid = 1L), get_individual(test_pop, pid = 5L))) expect_equal(-1L, meiotic_dist(get_individual(test_pop, pid = 37L), get_individual(test_pop, pid = 10L))) }) test_that("meiotic_dist_threshold works", { expect_equal(4L, meiotic_dist_threshold(get_individual(test_pop, pid = 1L), get_individual(test_pop, pid = 10L), threshold = 5)) expect_equal(4L, meiotic_dist_threshold(get_individual(test_pop, pid = 1L), get_individual(test_pop, pid = 10L), threshold = 4)) expect_equal(-1L, meiotic_dist_threshold(get_individual(test_pop, pid = 1L), get_individual(test_pop, pid = 10L), threshold = 3)) expect_equal(-1L, meiotic_dist_threshold(get_individual(test_pop, pid = 1L), get_individual(test_pop, pid = 10L), threshold = 2)) expect_equal(6L, meiotic_dist_threshold(get_individual(test_pop, pid = 1L), get_individual(test_pop, pid = 5L), threshold = 7)) expect_equal(6L, meiotic_dist_threshold(get_individual(test_pop, pid = 1L), get_individual(test_pop, pid = 5L), threshold = 6)) expect_equal(-1L, meiotic_dist_threshold(get_individual(test_pop, pid = 1L), get_individual(test_pop, pid = 5L), threshold = 5)) expect_equal(-1L, meiotic_dist_threshold(get_individual(test_pop, pid = 37L), get_individual(test_pop, pid = 10L), threshold = 1)) expect_equal(-1L, meiotic_dist_threshold(get_individual(test_pop, pid = 37L), get_individual(test_pop, pid = 10L), threshold = 0)) }) test_that("radius works 1", { if (FALSE) { plot(peds) } ans <- as.data.frame(meiotic_radius(get_individual(test_pop, pid = 1L), radius = 1)) expect_equal(nrow(ans), 2) expect_equal(subset(ans, pid == 1, dist, drop = TRUE), 0L) expect_equal(subset(ans, pid == 6, dist, drop = TRUE), 1L) ans <- as.data.frame(meiotic_radius(get_individual(test_pop, pid = 1L), radius = 2)) expect_equal(nrow(ans), 3) expect_equal(subset(ans, pid == 1, dist, drop = TRUE), 0L) expect_equal(subset(ans, pid == 6, dist, drop = TRUE), 1L) expect_equal(subset(ans, pid == 9, dist, drop = TRUE), 2L) ans <- as.data.frame(meiotic_radius(get_individual(test_pop, pid = 1L), radius = 3)) expect_equal(nrow(ans), 5) expect_equal(subset(ans, pid == 1, dist, drop = TRUE), 0L) expect_equal(subset(ans, pid == 6, dist, drop = TRUE), 1L) expect_equal(subset(ans, pid == 9, dist, drop = TRUE), 2L) expect_equal(subset(ans, pid == 11, dist, drop = TRUE), 3L) expect_equal(subset(ans, pid == 7, dist, drop = TRUE), 3L) ans <- as.data.frame(meiotic_radius(get_individual(test_pop, pid = 1L), radius = 4)) expect_equal(nrow(ans), 8) expect_equal(subset(ans, pid == 1, dist, drop = TRUE), 0L) expect_equal(subset(ans, pid == 6, dist, drop = TRUE), 1L) expect_equal(subset(ans, pid == 9, dist, drop = TRUE), 2L) expect_equal(subset(ans, pid == 11, dist, drop = TRUE), 3L) expect_equal(subset(ans, pid == 7, dist, drop = TRUE), 3L) expect_equal(subset(ans, pid == 2, dist, drop = TRUE), 4L) expect_equal(subset(ans, pid == 3, dist, drop = TRUE), 4L) expect_equal(subset(ans, pid == 10, dist, drop = TRUE), 4L) ans <- as.data.frame(meiotic_radius(get_individual(test_pop, pid = 1L), radius = 5)) expect_equal(nrow(ans), 9) expect_equal(subset(ans, pid == 1, dist, drop = TRUE), 0L) expect_equal(subset(ans, pid == 6, dist, drop = TRUE), 1L) expect_equal(subset(ans, pid == 9, dist, drop = TRUE), 2L) expect_equal(subset(ans, pid == 11, dist, drop = TRUE), 3L) expect_equal(subset(ans, pid == 7, dist, drop = TRUE), 3L) expect_equal(subset(ans, pid == 2, dist, drop = TRUE), 4L) expect_equal(subset(ans, pid == 3, dist, drop = TRUE), 4L) expect_equal(subset(ans, pid == 10, dist, drop = TRUE), 4L) expect_equal(subset(ans, pid == 8, dist, drop = TRUE), 5L) ##### ans <- as.data.frame(meiotic_radius(get_individual(test_pop, pid = 37L), radius = 1)) expect_equal(nrow(ans), 5) expect_equal(subset(ans, pid == 39, dist, drop = TRUE), 1L) expect_equal(subset(ans, pid == 33, dist, drop = TRUE), 1L) expect_equal(subset(ans, pid == 34, dist, drop = TRUE), 1L) expect_equal(subset(ans, pid == 35, dist, drop = TRUE), 1L) ans <- as.data.frame(meiotic_radius(get_individual(test_pop, pid = 37L), radius = 2)) expect_equal(nrow(ans), 7) expect_equal(subset(ans, pid == 39, dist, drop = TRUE), 1L) expect_equal(subset(ans, pid == 33, dist, drop = TRUE), 1L) expect_equal(subset(ans, pid == 34, dist, drop = TRUE), 1L) expect_equal(subset(ans, pid == 35, dist, drop = TRUE), 1L) expect_equal(subset(ans, pid == 36, dist, drop = TRUE), 2L) expect_equal(subset(ans, pid == 38, dist, drop = TRUE), 2L) ans <- as.data.frame(meiotic_radius(get_individual(test_pop, pid = 37L), radius = 3)) expect_equal(nrow(ans), 9) expect_equal(subset(ans, pid == 39, dist, drop = TRUE), 1L) expect_equal(subset(ans, pid == 33, dist, drop = TRUE), 1L) expect_equal(subset(ans, pid == 34, dist, drop = TRUE), 1L) expect_equal(subset(ans, pid == 35, dist, drop = TRUE), 1L) expect_equal(subset(ans, pid == 36, dist, drop = TRUE), 2L) expect_equal(subset(ans, pid == 38, dist, drop = TRUE), 2L) expect_equal(subset(ans, pid == 31, dist, drop = TRUE), 3L) expect_equal(subset(ans, pid == 32, dist, drop = TRUE), 3L) }) test_that("radius works 2", { if (FALSE) { plot(peds) ans <- meiotic_radius(get_individual(test_pop, pid = 1L), radius = 1) ans ans <- meiotic_radius(get_individual(test_pop, pid = 1L), radius = 2) ans ans <- meiotic_radius(get_individual(test_pop, pid = 1L), radius = 3) ans ans <- meiotic_radius(get_individual(test_pop, pid = 1L), radius = 4) ans } verify_distances <- function(pid_PoI, radius) { ans <- meiotic_radius(get_individual(test_pop, pid = pid_PoI), radius = radius) ans expect_true(isTRUE(all(ans[, 2] <= radius))) for (i in seq_len(nrow(ans))) { expect_equal(as.integer(ans[i, 2L]), meiotic_dist(get_individual(test_pop, pid = pid_PoI), get_individual(test_pop, pid = as.integer(ans[i, 1L]))), label = paste0("pid_PoI = ", pid_PoI, " / radius = ", radius, " / pid ", as.integer(ans[i, 1L]))) } } for (radius in seq_len(10)) { verify_distances(pid_PoI = 1L, radius = radius) } for (radius in seq_len(10)) { verify_distances(pid_PoI = 34L, radius = radius) } }) test_that("radius works 3", { set.seed(20230212) radius_pop <- sample_geneology(population_size = 1e2, generations = 30, generations_full = 3, generations_return = 3, progress = FALSE) radius_peds <- build_pedigrees(radius_pop$population, progress = FALSE) if (FALSE) { plot(radius_peds) } radius_verify_distances <- function(pid_PoI, radius) { #cat("pid_PoI = ", pid_PoI, " / radius = ", radius, "\n", sep = "") ans <- meiotic_radius(get_individual(radius_pop$population, pid = pid_PoI), radius = radius) ans expect_true(isTRUE(all(ans[, 2] <= radius))) for (i in seq_len(nrow(ans))) { expect_equal(as.integer(ans[i, 2L]), meiotic_dist(get_individual(radius_pop$population, pid = pid_PoI), get_individual(radius_pop$population, pid = as.integer(ans[i, 1L]))), label = paste0("pid_PoI = ", pid_PoI, " / radius = ", radius, " / pid ", as.integer(ans[i, 1L]))) } } radius_pop_live_pids <- lapply(radius_pop$individuals_generations, get_pid) pid_PoIs <- sample(x = radius_pop_live_pids, size = 5) for (pid_PoI in pid_PoIs) { for (radius in seq_len(15)) { radius_verify_distances(pid_PoI = pid_PoI, radius = radius) } } }) LOCI <- 5L pedigrees_all_populate_haplotypes(peds, loci = LOCI, mutation_rates = rep(0, LOCI), progress = FALSE) test_that("haplotype_matches_individuals works", { expect_equal(length(indvs), length(haplotype_matches_individuals(indvs, rep(0L, LOCI)))) expect_equal(lapply(indvs, get_pid), lapply(haplotype_matches_individuals(indvs, rep(0L, LOCI)), get_pid)) }) ################################################ indvs <- get_individuals(test_pop) set.seed(1) pedigrees_all_populate_haplotypes(peds, loci = LOCI, mutation_rates = rep(0.1, LOCI), progress = FALSE) test_that("count_haplotype_occurrences_individuals works", { expect_true(0L < count_haplotype_occurrences_individuals(indvs, rep(0L, LOCI))) }) test_that("get_matching_pids_from_hashmap works for (0, 0, ..., 0)", { hashmap <- build_haplotype_hashmap(indvs, progress = FALSE) pids <- get_matching_pids_from_hashmap(hashmap, rep(0L, LOCI)) expect_equal(length(pids), count_haplotype_occurrences_individuals(indvs, rep(0L, LOCI))) delete_haplotypeids_hashmap(hashmap) }) test_that("get_matching_pids_from_hashmap works for all", { set.seed(1) pedigrees_all_populate_haplotypes(peds, loci = LOCI, mutation_rates = rep(0.5, LOCI), progress = FALSE) all_hs <- get_haplotypes_individuals(indvs) all_hs <- all_hs[!duplicated(all_hs), ] hashmap <- build_haplotype_hashmap(indvs, progress = FALSE) for (i in seq_len(nrow(all_hs))) { h <- all_hs[i, ] pids <- get_matching_pids_from_hashmap(hashmap, h) expect_equal(length(pids), count_haplotype_occurrences_individuals(indvs, h), info = paste0("i = ", i, "; h = (", paste0(h, collapse = ", "), ")")) } delete_haplotypeids_hashmap(hashmap) })