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 r_cibd.

ibd

Integer vector with IBD states in segments if x is not a result from r_cibd. Taking values 0, 1, 2 for states = "ibd" or states = "kappa", 1, ..., 9 for states="identity" and 1, ..., 15 for states = "detailed".

pedigree

Pedigree in pedtools::ped form.

ids

Ids for which IBD is observed. Default is pedtools::leaves(pedigree).

states

One of "ibd" (default), "kappa", "identity" or "detailed".

log10

Should the log10 probability density be returned? Default is FALSE.

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]