### Trivial examples ###
x = nuclearPed(father = "A", mother = "B", children = "C")
# Random, distinct
patt1 = gip(x, list(c("A", "B"), c("A", "C")))
patt1
# Random, non-distinct
patt2 = gip(x, list(c("A", "B"), c("A", "C")), distinct = FALSE)
patt2
# Fully deterministic, distinct
patt3 = gip(x, list(c(p="A", p="C"), c(m="C")))
patt3
# Partially deterministic, non-distinct`
patt4 = gip(x, list(c("A", p="C"), c("B", m="C")), distinct = FALSE)
patt4
# Partially deterministic, single block
patt5 = gip(x, list(c(p="A", "A", "A")))
patt5
stopifnot(
gKinship(x, patt1) == 0, # (since A and B are unrelated)
gKinship(x, patt2) == 0, # (same as previous)
gKinship(x, patt3) == 0.5, # (only uncertainty is which allele A gave to C)
gKinship(x, patt4) == 0.25, # (distinct irrelevant)
gKinship(x, patt5) == 0.25 # (both random must hit the paternal)
)
### Kappa coefficients via generalised kinship ###
# NB: Much less efficient than `kappaIBD()`; only for validation
kappa_from_gk = function(x, ids, method = "WL") {
fa1 = father(x, ids[1])
fa2 = father(x, ids[2])
mo1 = mother(x, ids[1])
mo2 = mother(x, ids[2])
GK = function(...) gKinship(x, list(...), method = method)
k0 = GK(fa1, fa2, mo1, mo2)
k1 = GK(c(fa1, fa2), mo1, mo2) + GK(c(fa1, mo2), fa2, mo1) +
GK(c(mo1, fa2), fa1, mo2) + GK(c(mo1, mo2), fa1, fa2)
k2 = GK(c(fa1, fa2), c(mo1, mo2)) + GK(c(fa1, mo2), c(mo1, fa2))
c(k0, k1, k2)
}
y1 = nuclearPed(2); ids = 3:4
stopifnot(kappa_from_gk(y1, ids) == kappaIBD(y1, ids))
y2 = quadHalfFirstCousins(); ids = 9:10
stopifnot(kappa_from_gk(y2, ids) == kappaIBD(y2, ids))
### Detailed outputs and debugging ###
x = fullSibMating(1)
# Probability of sampling IBD alleles from 1, 5 and 6
p1 = gip(x, list(c(1,5,6)))
p1
gKinship(x, p1, method = "K", verbose = TRUE, debug = TRUE)
gKinship(x, p1, method = "WL", verbose = TRUE, debug = TRUE)
# Probability that paternal of 5 is IBD with maternal of 6
p2 = gip(x, list(c(p=5, m=6)))
p2
gKinship(x, p2, method = "LS", verbose = TRUE, debug = TRUE)
gKinship(x, p2, method = "GC", verbose = TRUE, debug = TRUE)
# Probability that paternal of 5 is *not* IBD with maternal of 6
p3 = gip(x, list(c(p=5), c(m=6)), distinct = TRUE)
p3
gKinship(x, p3, method = "LS", verbose = TRUE, debug = TRUE)
Run the code above in your browser using DataLab