d_cibd {ibdsegments} | R Documentation |
Probability density of continuously observed IBD segment lengths for pedigree members
Description
The d_cibd
function computes the probability density of observed
IBD segments on a chromosome.
Usage
d_cibd(
x,
ibd,
pedigree,
ids = pedtools::leaves(pedigree),
states = "ibd",
log10 = FALSE
)
Arguments
x |
Numeric vector with lengths of segments (centiMorgan) or result from |
ibd |
Integer vector with IBD states in segments if |
pedigree |
Pedigree in |
ids |
Ids for which IBD is observed. Default is |
states |
One of |
log10 |
Should the log10 probability density be returned? Default is |
Value
Numeric
Examples
ped_fs <- pedtools::nuclearPed(nch = 2)
# probability that full siblings are double IBD (kappa2)
d_cibd(x = 0, ibd = 2, ped_fs)
# full siblings are double IBD and remain so for 100 cM
d_cibd(x = 100, ibd = 2, ped_fs)
# full siblings are double IBD for 50 cM,
# then single IBD for 50 cM
d_cibd(x = c(50, 50), ibd = c(2, 1), ped_fs)
# full siblings are double IBD, remain so for 100 cM
# and no longer
d_cibd(x = c(100, 0), ibd = c(2, 1), ped_fs)
## probability density of IBD segment length for first cousins on an infinite chromosome
ped_fc <- pedtools::cousinPed()
# first compute the probability of IBD
k1_fc <- d_ibd(ibd = 1, ped_fc)
# density of segment length
f <- Vectorize(function(l) d_cibd(x = c(l,0), ibd = c(1, 0), ped_fc) / k1_fc)
curve(f, 0, 300)
# f is a probability density (integrates to 1)
integrate(f, 0, Inf)
# for full siblings, how does the chance of remaining double IBD
# depend on the segment length?
cM <- seq(from = 0, to = 100, length = 200)
pr_2ibd <- sapply(cM, d_cibd, 2, ped_fs) / d_ibd(2, ped_fs)
plot(cM, pr_2ibd, type="l")
# since there are four meioses, the sojourn time in this IBD state
# follows an exponential distribution with rate 0.04
stopifnot(all.equal(pr_2ibd, pexp(cM, rate = 0.04, lower.tail = FALSE)))
## Using the output from r_cibd directly to compute evidential value
## of IBD segments on chromosome for distinguishing first and second cousins
ped_c1 <- pedtools::cousinPed()
ped_c2 <- pedtools::cousinPed(degree = 2)
r_c1 <- r_cibd(n = 1e2, pedigree = ped_c1)
lr <- d_cibd(r_c1, pedigree = ped_c1) / d_cibd(r_c1, pedigree = ped_c2)
hist(log10(lr))
[Package ibdsegments version 1.0.1 Index]