## Example from Donnelly (1983)
# (historic) chromosome lengths (cM) used in Donnelly (1983)
L <- 33 * c(9.12, 8.53, 7.16, 6.59, 6.15, 5.87, 5.31, 4.92, 4.81, 4.71, 4.60,
4.47, 3.56, 3.60, 3.40, 3.20, 3.12, 2.72, 2.48, 2.27, 1.77, 1.64)
# Reproduce the "Offspring" column in Table 1
number_of_children <- 1:16
pr_survival <- pr_all_genes_survive(number_of_children, chromosome_length = L)
# Plot the results (compare to Fig. 10)
plot(number_of_children, pr_survival)
# Add the Poisson approximation
k <- 2:17 - 1
lines(k+1, exp(-k*33 / (2^k)), lty=2)
## Example using general purpose methods
# The probability of passing down all genes to k offspring is equal to the
# probability that the joint IBD state of k half siblings is 0 everywhere -
# there is no point where they all inherited the same DNA from the common
# parent.
# Define a function to create a pedigree of half siblings
ped_halfsibs <- function(number_of_half_sibs){
ped <- pedtools::nuclearPed(nch = 1)
for (k in 2:number_of_half_sibs) {
ped <- pedtools::addChild(ped, parents = c(1), verbose = FALSE)
}
ped
}
# Compute the probability that a chromosome of length 100 cm survives
# if the next generation consists of 5 children
p1 <- pr_all_genes_survive(number_of_children = 5L, chromosome_length = 100)
p2 <- d_cibd(x = 100, ibd = 0, pedigree = ped_halfsibs(5))
p1
p2
stopifnot(all.equal(p1, p2))
# If you have five children, which fraction of chromosome 1
# is passed down to the next generation? Assume a length of 268 cM
# Pr. of passing down the complete grand-maternal and grand-paternal parts
# is 42%
pr_all_genes_survive(number_of_children = 5L, chromosome_length = 268)
# The distribution of the fraction that is passed down has a point
# mass of 0.42 at 100% and has a continuous density with weight 0.58
dist <- total_ibd_dist(ped_halfsibs(5), ibd_state = 0,
chromosome_length = 267, fraction = TRUE)
dist
plot(dist)
Run the code above in your browser using DataLab