test-mixtures.R
In malan: MAle Lineage ANalysis

context("Mixtures")

set.seed(1)
sim_res_growth <- sample_geneology_varying_size(population_sizes = rep(1e3, 200),
                                                generations_full = 3,
                                                generations_return = 3,
                                                enable_gamma_variance_extension = TRUE,
                                                gamma_parameter_shape = 5,
                                                gamma_parameter_scale = 1/5,
                                                progress = FALSE)

pedigrees <- build_pedigrees(sim_res_growth$population, progress = FALSE)

mutrts <- rep(0.001, 20)
pedigrees_all_populate_haplotypes(pedigrees = pedigrees, 
                                  loci = length(mutrts), 
                                  mutation_rates = mutrts, progress = FALSE)
live_individuals <- sim_res_growth$individuals_generations

U_indices <- sample.int(n = length(live_individuals), size = 2, replace = FALSE)
U1 <- live_individuals[[U_indices[1]]]
U2 <- live_individuals[[U_indices[2]]]
H1 <- get_haplotype(U1)
H2 <- get_haplotype(U2)

# Max 100 redraws
for (i in 1:100) {
  if (any(H1 != H2)) {
    break
  }
  
  # If H1 == H2, try again:
  U_indices <- sample.int(n = length(live_individuals), size = 2, replace = FALSE)
  U1 <- live_individuals[[U_indices[1]]]
  U2 <- live_individuals[[U_indices[2]]]
  H1 <- get_haplotype(U1)
  H2 <- get_haplotype(U2)
}

test_that("contributors different", {
  expect_true(any(H1 != H2))
})


mixres <- mixture_info_by_individuals_2pers(live_individuals, U1, U2)
others_haps <- get_haplotypes_pids(sim_res_growth$population, mixres$pids_others_included)

others_indv <- lapply(mixres$pids_others_included, function(pid) {
  get_individual(sim_res_growth$population, pid)
})
others_haps_2 <- get_haplotypes_individuals(individuals = others_indv)

test_that("mixture_info_by_individuals_2pers works", {
  expect_equal(mixres$donor1_profile, H1)
  expect_equal(mixres$donor2_profile, H2)
  expect_equal(mixres$loci_not_matching, sum(H1 != H2))
  expect_equal(length(mixres$pids_matching_donor1),
               count_haplotype_occurrences_individuals(individuals = live_individuals, haplotype = H1))
  expect_equal(length(mixres$pids_matching_donor2),
               count_haplotype_occurrences_individuals(individuals = live_individuals, haplotype = H2))
  expect_equal(others_haps, others_haps_2)
})

others_haps_unique <- others_haps[!duplicated(others_haps), ]
others_haps_counts <- unlist(lapply(seq_len(nrow(others_haps_unique)), function(hap_i) {
  count_haplotype_occurrences_individuals(individuals = live_individuals, 
                                          haplotype = others_haps_unique[hap_i, ])
}))

test_that("mixture others included haplotypes works", {
  expect_equal(sum(others_haps_counts), length(mixres$pids_others_included))
})

##############################################################################

set.seed(1)
sim_res_growth <- sample_geneology_varying_size(population_sizes = rep(1e3, 200),
                                                generations_full = 3,
                                                generations_return = 3,
                                                enable_gamma_variance_extension = FALSE,
                                                progress = FALSE)

pedigrees <- build_pedigrees(sim_res_growth$population, progress = FALSE)

mutrts <- rep(0.1, 2)
pedigrees_all_populate_haplotypes_ladder_bounded(pedigrees = pedigrees, 
                                                 mutation_rates = mutrts, 
                                                 ladder_min = c(1, 1),
                                                 ladder_max = c(2, 3),
                                                 get_founder_haplotype = function() c(1, 1),
                                                 progress = FALSE)

live_individuals <- sim_res_growth$individuals_generations
indvs_haps <- get_haplotypes_individuals(live_individuals)

U1 <- live_individuals[[ which(apply(indvs_haps, 1, function(x) identical(x, c(1L, 1L))))[1L] ]]
U2 <- live_individuals[[ which(apply(indvs_haps, 1, function(x) identical(x, c(2L, 2L))))[1L] ]]
U3 <- live_individuals[[ which(apply(indvs_haps, 1, function(x) identical(x, c(2L, 3L))))[1L] ]]


###############################################################################################
# 2 persons
###############################################################################################
# mix = U1 (1, 1) + U2 (2, 2) = {1, 2} x {1, 2}
mixres2 <- mixture_info_by_individuals_2pers(live_individuals, U1, U2)

# Haps included in mixture
haps2_unique <- list(c(1L, 1L),  # <- U1 = SUSPECT / KNOWN
                     c(1L, 2L), 
                     c(2L, 1L), 
                     c(2L, 2L))
haps2_unique_counts <- c(11, 12, 21, 22) # counts refering to types
res2 <- analyse_mixture_result(mixres2, haps2_unique, haps2_unique_counts)

test_that("mixture 2 persons", {
  # Hp: only (2, 2) with count 22 must be explained:
  expect_equal(res2$terms_Hp, list(22))
  
  # Hd: mixture either (1, 1)+(2, 2) or (1, 2)+(2, 1) [(2, 1) + (1, 2) is same but considering order]
  # Two random males would give:
  #   expect_equal(res2$terms_Hd, list(c(11, 22), c(12, 21)))
  # But now we say two random different from the PoI, so substract by 1 for PoI's profile (11):
  expect_equal(res2$terms_Hd, list(c(11-1, 22), c(12, 21)))
})


###############################################################################################
# 3 persons
###############################################################################################
# mix = U1 (1, 1) + U2 (2, 2) + U3 (2, 3) = {1, 2} x {1, 2, 3}
mixres3 <- mixture_info_by_individuals_3pers(live_individuals, U1, U2, U3)

if (FALSE) {
  haps_in_mixture_3pers <- get_haplotypes_pids(sim_res_growth$population, mixres3$pids_included_in_mixture)
  haps_in_mixture_list_3pers <- lapply(seq_len(nrow(haps_in_mixture_3pers)), function(i) haps_in_mixture_3pers[i, , drop = FALSE])
  unique_haps_in_mixture_3pers <- haps_in_mixture_list_3pers[!duplicated(haps_in_mixture_3pers)]
  hapcount_pop_3gen_3pers <- unlist(lapply(seq_along(unique_haps_in_mixture_3pers), function(hap_i) {
    count_haplotype_occurrences_individuals(individuals = live_individuals, 
                                            haplotype = unique_haps_in_mixture_3pers[[hap_i]])
  }))
  do.call(rbind, unique_haps_in_mixture_3pers)
}

# Haps included in mixture, pop made such that this is {1, 2} x {1, 2, 3}
haps3_unique <- list(c(1L, 1L),  # <- U1 = SUSPECT / KNOWN
                     c(1L, 2L), 
                     c(1L, 3L),
                     c(2L, 1L), 
                     c(2L, 2L),
                     c(2L, 3L))
haps3_unique_counts <- c(11, 12, 13, 21, 22, 23) # counts refering to types
res3 <- analyse_mixture_result(mixres3, haps3_unique, haps3_unique_counts)

test_that("mixture 3 persons", {
  # Hp: allele 2 at locus 1 and alleles {2, 3} at locus 2 must be explained.
  #     thus unknowns are bound to have 2/3 at locus 2, but first locus can change, 
  #     at least one must have 2 at locus 1, but both can potentially:
  #     (L1, L2)
  #     (1, 2) + (2, 3)
  #     (2, 2) + (1, 3)
  #     (2, 2) + (2, 3)
  expect_equal(res3$terms_Hp, list(
    c(12, 23), 
    c(13, 22),
    c(22, 23)
    ))
  
  # Hd: 
  #     (L1, L2)
  #     (1, 1) + (1, 2) + (2, 3)
  #     (1, 1) + (2, 2) + (2, 3)
  #     (1, 1) + (2, 2) + (1, 3) -> [enumeration order] -> (1, 1) + (1, 3) + (2, 2)
  # 
  #     (1, 2) + (1, 3) + (2, 1)
  #     (1, 2) + (2, 1) + (2, 3)
  # 
  #     (1, 3) + (2, 1) + (2, 2)
  
  # Three random males would give:
  #   expect_equal(res3$terms_Hd, list(
  #     c(11, 12, 23), 
  #     c(11, 13, 22), # (1, 1) + (2, 2) + (1, 3) -> [enumeration order] -> (1, 1) + (1, 3) + (2, 2)
  #     c(11, 22, 23),
  #   
  #     c(12, 13, 21),
  #     c(12, 21, 23),
  #   
  #     c(13, 21, 22)
  #   ))
  #   expect_equal(res2$terms_Hd, list(c(11, 22), c(12, 21)))
  # But now we say three random different from the PoI, so substract by 1 for PoI's profile (11):
  expect_equal(res3$terms_Hd, list(
    c(11-1, 12, 23), 
    c(11-1, 13, 22), # (1, 1) + (2, 2) + (1, 3) -> [enumeration order] -> (1, 1) + (1, 3) + (2, 2)
    c(11-1, 22, 23),
    
    c(12, 13, 21),
    c(12, 21, 23),
    
    c(13, 21, 22)
  ))
})