quickEP = function(claim, true, ids, afreq = NULL, ...) { als = if(is.null(afreq)) NULL else seq_along(afreq) exclusionPower(claim, true, ids, alleles = als, afreq = afreq, plot = F, verbose = F, ...)$EPtotal } test_that("EP works in empty paternity case", { claim = nuclearPed(1) true = singletons(c(1,3)) ids = c(1, 3) afr = c(.1, .9) # Autosomal ep_aut = quickEP(claim, true, ids, afr) expect_equal(ep_aut, 2 * afr[1]^2 * afr[2]^2) # X: father/son (no exclusion) expect_equal(quickEP(claim, true, ids, afr, Xchrom = T), 0) # X: father/daughter claim2 = nuclearPed(1, sex = 2) true2 = singletons(c(1,3), sex = 1:2) ep_X = quickEP(claim2, true2, ids, afr, Xchrom = T) expect_equal(ep_X, afr[1] * afr[2]) # p^2*q + p*q^2 }) test_that("EP works in empty pat-case with added singletons", { claim = list(nuclearPed(1), singleton(4)) true = singletons(c(1,3,4)) ids = c(1, 3, 4) afr = c(.1, .9) # Autosomal expect_equal(quickEP(claim, true, ids, afr), 2 * afr[1]^2 * afr[2]^2) # X expect_equal(quickEP(claim, true, ids, afr, Xchrom = T), 0) # X: father/daughter claim2 = list(nuclearPed(1, sex = 2), singleton(4)) true2 = singletons(c(1,3,4), sex = c(1,2,1)) expect_equal(quickEP(claim2, true2, ids, afr, Xchrom = T), afr[1] * afr[2]) }) test_that("EP works in paternity case with child typed", { claim = nuclearPed(1, sex = 2) true = singletons(c(1,3)) afr = c(.5, .3, .2) m = mX = marker(claim, `3` = 1, alleles = 1:3, afreq = afr) chrom(mX) = 23L claim = setMarkers(claim, list(m, mX)) # Autosomal ep_aut = quickEP(claim, true, ids = 1, markers = 1) expect_equal(ep_aut, sum(afr[-1])^2) # X ep_X = quickEP(claim, true, ids = 1, markers = 2) expect_equal(ep_X, sum(afr[-1])) }) test_that("EP works in paternity case with parents typed", { claim = nuclearPed(1, sex = 2) true = singletons(c(1,3), sex = 1:2) afr = c(.5, .3, .2) m = mX = marker(claim, `1` = 1, `2` = 2, alleles = 1:3, afreq = afr) chrom(mX) = 23L claim = setMarkers(claim, list(m, mX)) # Autosomal ep_aut = quickEP(claim, true, ids = 3, markers = 1) expect_equal(ep_aut, 1 - 2*afr[1]*afr[2]) # X ep_X = quickEP(claim, true, ids = 3, markers = 2) expect_equal(ep_X, 1 - 2*afr[1]*afr[2]) })