Title: | Identity by Descent Probability in Pedigrees |
Version: | 1.0.1 |
Description: | Identity by Descent (IBD) distributions in pedigrees. A Hidden Markov Model is used to compute identity coefficients, simulate IBD segments and to derive the distribution of total IBD sharing and segment count across chromosomes. The methods are applied in Kruijver (2025) <doi:10.3390/genes16050492>. The probability that the total IBD sharing is zero can be computed using the method of Donnelly (1983) <doi:10.1016/0040-5809(83)90004-7>. |
License: | GPL-2 | GPL-3 [expanded from: GPL (≥ 2)] |
Encoding: | UTF-8 |
RoxygenNote: | 7.3.2 |
LinkingTo: | Rcpp |
Imports: | Rcpp, pedtools, expm |
Suggests: | testthat (≥ 3.0.0), ribd, ggplot2 |
Config/testthat/edition: | 3 |
NeedsCompilation: | yes |
Packaged: | 2025-04-29 16:50:03 UTC; maart |
Author: | Maarten Kruijver |
Maintainer: | Maarten Kruijver <maarten.kruijver@esr.cri.nz> |
Repository: | CRAN |
Date/Publication: | 2025-05-02 09:30:05 UTC |
ibdsegments: Identity by Descent Probability in Pedigrees
Description
Identity by Descent (IBD) distributions in pedigrees. A Hidden Markov Model is used to compute identity coefficients, simulate IBD segments and to derive the distribution of total IBD sharing and segment count across chromosomes. The methods are applied in Kruijver (2025) doi:10.3390/genes16050492. The probability that the total IBD sharing is zero can be computed using the method of Donnelly (1983) doi:10.1016/0040-5809(83)90004-7.
Author(s)
Maintainer: Maarten Kruijver maarten.kruijver@esr.cri.nz (ORCID)
Expectation for IBD Distribution Objects
Description
Computes the expectation (mean) of a total IBD or segment count distribution.
Usage
E(x, ...)
Arguments
x |
An object of class |
... |
Additional arguments (not used). |
Value
A numeric value.
Convolve IBD distributions to obtain the distribution of the sum
Description
Chromosome-wise IBD distributions obtained from the total_ibd_dist
function can be convoluted manually using the convolve_total_ibd_dists
function. This allows finer control of the procedure by controlling
the number of gridpoints used in the FFT and the threshold for
point masses to be retained.
Usage
convolve_total_ibd_dists(
...,
point_mass_eps = 1e-09,
number_of_gridpoints_exponent = 12
)
Arguments
... |
ibd dists. |
point_mass_eps |
Point masses smaller than this value are discarded. |
number_of_gridpoints_exponent |
Default is 12. |
Details
The convolve_total_ibd_dists
implements a convolution procedure based on a self
contained variant of the recipe implemented in the distr
package adapted
to this use case. In particular, the IBD distribution for one
chromosome is often a mixture of two point masses and continuous with finite
support otherwise. Convolution increases the number of point masses and
decreases their probability mass. This function assumes that small point
masses are not of specific interest and will discard the point masses when
the probability mass is smaller than point_mass_eps
with a default value
of 1e-9
. For typical pedigree relationships this means that the IBD
distribution of more than a few chromosomes is treated as a continuous
distribution admitting a density function.
The number_of_gridpoints_exponent
controls the density of gridpoints in the
FFT. By default this is 12
which means that 2^12=4096
gridpoints are
used. Increasing this parameter improves accuracy at the cost of increased
runtime and memory use.
The point masses are of particular interest in some applications. For
example, the probability that no part of the autosomal genome is inherited
may be of interest. In that case, the point masses should not be discarded
and the d_cibd
function may be used. See the example below, for further
details.
Value
Object of class total_ibd_dist
Examples
## Convolution of IBD distributions for half siblings at chromosome 1 and 2
# define half sib pedigree
ped_hs <- pedtools::halfSibPed()
# obtain chromosome-wise distributions
d1 <- total_ibd_dist(pedigree = ped_hs, chromosome_length = 267.77)
d2 <- total_ibd_dist(pedigree = ped_hs, chromosome_length = 251.73)
convolve_total_ibd_dists(d1, d2) # 4 point masses
## Accuracy of convolution depends on number of gridpoints
L <- c(267.77, 251.73, 218.31, 202.89, 197.08, 186.02, 178.4, 161.54,
157.35, 169.28, 154.5, 165.49, 127.23, 116, 117.32, 126.59, 129.53,
116.52, 106.35, 107.76, 62.88, 70.84)
# obtain chromosome-wise distributions for half siblings
hs <- total_ibd_dist(pedigree = ped_hs,
chromosome_length = L, convolve = FALSE)
# obtain moments of the sum by summing the moments..
mean_hat <- sum(sapply(hs, E))
sd_hat <- sqrt(sum(sapply(hs, var)))
# .. or by summing the distributions with varying numbers of gridpoints
k <- 6:10
sd_hat_by_k <- sapply(k, function(k) sd(convolve_total_ibd_dists(hs,
number_of_gridpoints_exponent = k)))
mean_hat_by_k <- sapply(k, function(k) E(convolve_total_ibd_dists(hs,
number_of_gridpoints_exponent = k)))
plot(k, mean_hat_by_k)
abline(h = mean_hat, lty = 2)
plot(k, sd_hat_by_k)
abline(h = sd_hat, lty = 2)
Probability function for IBD Distribution Objects
Description
Returns the probability function of a total IBD or segment count distribution.
Usage
d(x, ...)
Arguments
x |
An object of class |
... |
Additional arguments (not used). |
Value
A numeric value.
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))
Compute probability of IBD for pedigree members
Description
The d_ibd
function computes the likelihood of IBD
for one position or multiple linked markers on the same chromosome.
Usage
d_ibd(
ibd,
pedigree,
ids = pedtools::leaves(pedigree),
recombination_rate_by_locus = numeric(),
states = "ibd",
log10 = FALSE
)
Arguments
ibd |
Integer vector. Taking values 0, 1, 2 for |
pedigree |
Pedigree in |
ids |
Ids for which IBD is observed. Defaults to |
recombination_rate_by_locus |
Numeric vector with length one shorter than |
states |
One of |
log10 |
Should the log10 likelihood be returned? Default is |
Value
Numeric
Examples
# Compute kappa0, kappa1, kappa2 for full siblings
ped_fs <- pedtools::nuclearPed(nch = 2)
k0 <- d_ibd(ibd = 0, pedigree = ped_fs)
k1 <- d_ibd(ibd = 1, pedigree = ped_fs)
k2 <- d_ibd(ibd = 2, pedigree = ped_fs)
c(k0, k1, k2)
stopifnot(identical(c(k0, k1, k2), c(0.25, 0.5, 0.25)))
# Compute kappa00 for two tightly linked loci
d_ibd(c(0,0), pedigree = ped_fs,
recombination_rate_by_locus = c(0.01))
# or 100 tightly linked loci
d_ibd(rep(0, 100), pedigree = ped_fs,
recombination_rate_by_locus = c(rep(0.01, 99)))
# Jacquard's 9 condensed and 15 detailed identity coefficients
ped_fs_mating <- pedtools::fullSibMating(1)
sapply(1:9, d_ibd, pedigree = ped_fs_mating, states = "identity")
sapply(1:15, d_ibd, pedigree = ped_fs_mating, states = "detailed")
Inheritance pattern for inheritance vectors
Description
The inheritance_pattern
function determines the inheritance pattern for all
pedigree members by dropping the founder allele labels down the pedigree
according to the IBD vector v
.
Usage
inheritance_pattern(inheritance_space, v)
Arguments
inheritance_space |
Output of inheritance_space. |
v |
Integer. |
Value
Data frame describing inheritance pattern and IBD state for each v
.
Examples
ped_fs <- pedtools::nuclearPed(nch = 2)
i <- inheritance_space(ped_fs, ids = 3:4)
# show the inheritance pattern and IBD state for all canonical IBD vectors
inheritance_pattern(i, v = 0:3)
# without exploiting founder symmetry
i2 <- inheritance_space(ped_fs, ids = 3:4, exploit_symmetries = FALSE)
inheritance_pattern(i2, v = 0:15)
Inheritance space for pedigree
Description
The inheritance_space
function determines the space of IBD vectors for a pedigree.
This is mostly for internal use but may be interesting by itself.
Usage
inheritance_space(pedigree, ids, states = "ibd", exploit_symmetries = TRUE)
Arguments
pedigree |
Pedigree in |
ids |
Ids for which IBD is observed. Default is |
states |
One of |
exploit_symmetries |
Should symmetries be used to reduce to state space?
This can be set to |
Value
Object of class inheritance_space
.
Examples
# set up inheritance space for half sib pedigree
i <- inheritance_space(pedigree = pedtools::halfSibPed())
# since there are 2 non-founders, there are 2^4 IBD vectors
# but only 2 distinct states are considered because of symmetries
i
# pry into the internals to see individual pedigree transmissions
i$transmissions
Compute probabilities of minimal multi-id IBD patterns
Description
For two full siblings at one locus, it is well known that there are three distinct minimal IBD patterns
with probabilities 0.25, 0.5 and 0.25. The ribd::multiPersonIBD
function
generalises the computation of these patterns and their probabilities to more than
two ids The multi_ibd_patterns
function further generalises the computation
to patterns across multiple loci.
Usage
multi_ibd_patterns(
pedigree,
ids = pedtools::leaves(pedigree),
recombination_rate_by_locus = numeric()
)
Arguments
pedigree |
Pedigree in |
ids |
Ids for which IBD is observed. Defaults to |
recombination_rate_by_locus |
Optionally a numeric vector with recombination rates. |
Value
DataFrame
Examples
# Compute IBD patterns for two full siblings...
multi_ibd_patterns(pedtools::nuclearPed(nch = 2))
# ... and the generalisation to three siblings
multi_ibd_patterns(pedtools::nuclearPed(nch = 3))
# Two full siblings at two tightly linked loci
multi_ibd_patterns(pedtools::nuclearPed(nch = 2),
recombination_rate_by_locus = 0.01)
Probability of no IBD across one or more chromosomes
Description
Donnelly (1983) studies the probability that relatives of certain types have no segments of a chromosome in common and provides expressions that can be efficiently computed.
Usage
pr_0_total_ibd(
relationship_type = c("cousin", "grandparent", "halfsib", "uncle"),
degree,
removal,
removal1,
removal2,
chromosome_length,
log10 = FALSE
)
Arguments
relationship_type |
One of |
degree |
See details. |
removal |
See details. |
removal1 |
See details. |
removal2 |
See details. |
chromosome_length |
Default is 267.77 cM (an estimate of the length of chromosome 1). |
log10 |
Should the log10 probability be returned? Default is |
Details
Types of relationships supported are:
cousin
Use
degree = a
andremoval = b
fora
'th cousinsb
times removed with degree at least one. Default isdegree = 1
.halfsib
Use
removal = 0
(default) for half-siblings andremoval = a
for half-cousinsa
times removed.grandparent
degree = 1
(default) for grandparents, 2 for great-grandparents and so on.uncle
Use
degree = 0
(default) for uncle anddegree = d
for great^d uncle
Value
Numeric
References
Donnelly K. P. (1983). The probability that related individuals share some section of genome identical by descent. Theoretical population biology, 23(1), 34–63. https://doi.org/10.1016/0040-5809(83)90004-7
See Also
pr_all_genes_survive for the probability that all autosomal genes are passed on to the next generation (offspring column in Table 1 of Donnelly (1983))
Examples
## Cousin-type: third cousins on a chromosome of length 100 cM
degree <- 3
p0_3C <- pr_0_total_ibd("cousin", degree = 3, chromosome_length = 100)
p0_3C
# verify
p0_3C_manual <- d_cibd(x = 100, ibd = 0,
pedigree = pedtools::cousinPed(degree = 3))
p0_3C_manual
stopifnot(all.equal(p0_3C, p0_3C_manual))
## Half-sib type: half-cousins twice removed
p0_H1C_2R <- pr_0_total_ibd("halfsib",
degree = 1, removal = 2, chromosome_length = 100)
p0_H1C_2R
p0_H1C_2R_manual <- d_cibd(x = 100, ibd = 0,
pedigree = pedtools::halfCousinPed(removal = 2))
p0_H1C_2R_manual
stopifnot(all.equal(p0_H1C_2R, p0_H1C_2R_manual))
## Grandparent-type: great grandparents (degree = 2)
p0_GGP <- pr_0_total_ibd("grandparent", degree = 2, chromosome_length = 100)
p0_GGP
# GGP is a third generation ancestor so n = 3
p0_GGP_manual <- d_cibd(x = 100, ibd = 0,
pedigree = pedtools::linearPed(n = 3),
ids = c(1, pedtools::leaves(pedtools::linearPed(n = 3))))
p0_GGP_manual
stopifnot(all.equal(p0_GGP, p0_GGP_manual))
## Uncle-type: degree = 0 for uncle
p0_AV <- pr_0_total_ibd("uncle", chromosome_length = 100)
p0_AV
p0_AV_manual <- d_cibd(x = 100, ibd = 0, pedigree = pedtools::avuncularPed())
p0_AV_manual
stopifnot(all.equal(p0_AV, p0_AV_manual))
## Reproduce Table 1 of 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)
k <- 4:15
tab1 <- data.frame(k=k)
tab1$cousin <- pr_0_total_ibd(relationship_type = "cousin",
degree = rep(2:7, each = 2),
removal = rep(0:1, times = 6),
chromosome_length = L)
tab1$uncle <- pr_0_total_ibd(relationship_type = "uncle",
degree = k - 1, chromosome_length = L)
# Note the removal on one side only
tab1$halfsib <- pr_0_total_ibd(relationship_type = "halfsib",
removal1 = k - 1,
removal2 = rep(0, length(k)),
chromosome_length = L)
# Poisson approximation
tab1$`exp(-33 * k / 2^k)` <- exp(-33 * k / 2^k)
# Note that k corresponds to great^k grandparent,
# i.e. a (k+2)'th generation ancestor
# (not great^(k-1) and (k+1)'th generation ancestor as printed)
tab1$grandparent <- pr_0_total_ibd(relationship_type = "grandparent",
degree = k, chromosome_length = L)
tab1
Probability of passing down all autosomal genes to offspring
Description
Donnelly (1983) presents a way to efficiently calculate the probability of passing down all genes to the next generation. This includes both the paternally and maternally inherited genes, so at least two offspring are needed for this probability to be non-zero.
Usage
pr_all_genes_survive(
number_of_children,
chromosome_length = 267.66,
log10 = FALSE
)
Arguments
number_of_children |
How many offspring? Result is vectorised over this parameter. |
chromosome_length |
Default is 267.77 cM (an estimate of the length of
chromosome 1). |
log10 |
Should the log10 probability be returned? Default is |
Value
Numeric
References
Donnelly K. P. (1983). The probability that related individuals share some section of genome identical by descent. Theoretical population biology, 23(1), 34–63. https://doi.org/10.1016/0040-5809(83)90004-7
Examples
## 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)
Random generation for IBD on a continuous genome
Description
The r_cibd
function generates random Identity-by-Descent (IBD) segments
along a continuous genome, given a specified pedigree and observed
individuals.
Usage
r_cibd(
n,
pedigree,
ids = pedtools::leaves(pedigree),
states = "ibd",
ibd_state = 1L,
chromosome_length = 267.77
)
Arguments
n |
Number of observations |
pedigree |
Pedigree in |
ids |
Ids for which IBD is observed. Default is |
states |
One of |
ibd_state |
Default is 1. |
chromosome_length |
Default is 267.77 cM (an estimate of the length of chromosome 1). |
Value
A list containing:
samples |
Data frame of simulated IBD segments along the chromosome. |
stats |
Data frame with summary statistics per sample, including total IBD length and the segment count. |
Examples
## Basic example: IBD along one chromosome for half siblings
L <- 300
r_hs <- r_cibd(n = 1e4, pedigree = pedtools::halfSibPed(), chromosome_length = 300)
# half sibs alternate between IBD (state 1) and not IBD (state 0)
head(r_hs$samples)
# the total_length and number of segments per sample are also returned
head(r_hs$stats)
## Comparing half siblings and grandparent-grandchild
r_gp <- r_cibd(n = 1e4, pedigree = pedtools::linearPed(2), ids = c(1, 5),
chromosome_length = 300)
hist(r_gp$stats$total_length)
hist(r_hs$stats$total_length)
Standard Deviation for IBD Distribution Objects
Description
Computes the standard deviation of a total IBD or segment count distribution.
Usage
sd(x, ...)
Arguments
x |
An object of class |
... |
Additional arguments (not used). |
Value
A numeric value.
Compute probability distribution of IBD segment count
Description
The segment_count_dist
function computes the probability distribution
of the number of IBD segments
(i.e., the count of IBD intervals) between pairs of individuals
in a pedigree, for a given IBD state and chromosome length(s).
Usage
segment_count_dist(
pedigree,
ids = pedtools::leaves(pedigree),
states = "ibd",
ibd_state = 1L,
chromosome_length = 267.77,
convolve = TRUE
)
Arguments
pedigree |
Pedigree in |
ids |
Ids for which IBD is observed. Default is |
states |
One of |
ibd_state |
Default is 1. |
chromosome_length |
Default is 267.77 cM (an estimate of the length of chromosome 1). |
convolve |
Should the distribution of the sum (across chromosomes) be obtained? |
Details
This function is analogous to total_ibd_dist()
but focuses on
the number of IBD segments rather than the total IBD length.
Value
object of class segment_count_dist
See Also
total_ibd_dist()
for the distribution of IBD length.
Examples
ped_hs <- pedtools::halfSibPed()
dist <- segment_count_dist(ped_hs, chromosome_length = 100)
dist
plot(dist)
E(dist)
sd(dist)
r <- r_cibd(n = 1e4, pedigree = ped_hs, chromosome_length = 100)
mean(r$stats$segment_count)
sd(r$stats$segment_count)
Compute probability distribution of total IBD
Description
The total_ibd_dist
function computes the probability distribution
of the total IBD length (or fraction) over one or more autosomes.
Usage
total_ibd_dist(
pedigree,
ids = pedtools::leaves(pedigree),
fraction = FALSE,
states = "ibd",
ibd_state = 1L,
chromosome_length = 267.77,
convolve = TRUE,
...
)
Arguments
pedigree |
Pedigree in |
ids |
Ids for which IBD is observed. Default is |
fraction |
If TRUE, the distribution of the IBD fraction instead of length will be returned. Default is FALSE. |
states |
One of |
ibd_state |
Default is 1. |
chromosome_length |
Default is 267.77 cM (an estimate of the length of chromosome 1). |
convolve |
Should the distribution of the sum (across chromosomes) be obtained? |
... |
Additional parameters passed to |
Details
For many pedigree relationships, the probability distribution of the total IBD length over one chromosome is a mixture of two point masses (0 and chromosome length) and a continuous density.
If convolve=TRUE
(the default) and chromosome_length
has length
greater than one, the convolution of the distributions will be obtained
by FFT using the convolve_total_ibd_dists
function. Convolution will
typically produce a rapidly increasing number of point masses
with very small probabilities which are discarded if the
probability falls below a threshold of 1e-9
; see convolve_total_ibd_dists
for details and finer control.
Value
object of class total_ibd_dist
Examples
## Total IBD and fraction of IBD for a cousin relationship
ped_fc <- pedtools::cousinPed()
# total IBD length for 100 cM
dist_length <- total_ibd_dist(ped_fc, chromosome_length = 100)
dist_length
plot(dist_length)
# fraction IBD for 100 cM
dist_fraction <- total_ibd_dist(ped_fc, chromosome_length = 100,
fraction = TRUE)
dist_fraction
plot(dist_fraction)
# distribution of total length across three chromosomes (150, 200, 250 cM)
plot(total_ibd_dist(ped_fc, chromosome_length = c(150, 200, 250)))
# a quick approximation with reasonable accuracy (with just 256 gridpoints)
plot(total_ibd_dist(ped_fc, chromosome_length = c(150, 200, 250),
number_of_gridpoints_exponent = 8))
## Difference between IBD distributions between half-sibs, uncle-nephew
## and grandparent-grandchild relationships
# kappa1 is 1/2 for half sibs, uncle-nephew and grandparent-grandchild
# but the distribution of the fraction of a chromosome that is in this
# state differs between the relationships
# define pedigrees and verify kappa1
ped_gp <- pedtools::linearPed(n = 2)
ped_av <- pedtools::avuncularPed()
ped_hs <- pedtools::halfSibPed()
stopifnot(all.equal(1/2,
d_ibd(1, ped_av),
d_ibd(1, ped_hs),
d_ibd(1, ped_gp, ids = c(1,5))))
# Compute IBD distributions
d_av <- total_ibd_dist(ped_av)
d_hs <- total_ibd_dist(ped_hs)
d_gp <- total_ibd_dist(ped_gp, ids = c(1,5))
# the point masses are different
d_av
d_hs
d_gp
# plot the continuous densities
x0 <- seq(0, 267.77, length.out = 200)
df <- data.frame(cM = rep(x0, 3),
y = c(d(d_av)(x0), d(d_hs)(x0), d(d_gp)(x0)),
Relationship = rep(c("Avuncular", "Half-Sibling",
"Grandparent"), each = length(x0)))
require(ggplot2)
ggplot(df, aes(x = cM, y = y, color = Relationship)) + geom_line()
Compute moments of probability distribution of total IBD
Description
The total_ibd_dist_moments
function computes mean and variance of
the probability distribution of the total IBD length (or fraction)
over one autosome. The function uses double numerical integration.
Usage
total_ibd_dist_moments(
pedigree,
ids = pedtools::leaves(pedigree),
fraction = FALSE,
states = "ibd",
ibd_state = 1L,
chromosome_length = 267.77
)
Arguments
pedigree |
Pedigree in |
ids |
Ids for which IBD is observed. Default is |
fraction |
If TRUE, the distribution of the IBD fraction instead of length will be returned. Default is FALSE. |
states |
One of |
ibd_state |
Default is 1. |
chromosome_length |
Default is 267.77 cM (an estimate of the length of chromosome 1). |
Value
list
Examples
# Full Siblings are double IBD with 25% probability
# we may compute the expectation and variance of the total length of
# a chromosome that is double IBD for a chromosome of 100 cM (i.e. 1 Morgan)
m <- total_ibd_dist_moments(pedigree = pedtools::nuclearPed(nch = 2),
ibd_state = 2, chromosome_length = 100)
m
# compare to numerical integration of the full distribution
d_fs <- total_ibd_dist(pedigree = pedtools::nuclearPed(nch = 2),
ibd_state = 2, chromosome_length = 100)
m2 <- list(mean = E(d_fs), variance = var(d_fs), sd = sd(d_fs))
m2
stopifnot(all.equal(m, m2))
# Expectation and variance of _fraction_ of the genome that is
# double IBD between four full siblings
m4 <- total_ibd_dist_moments(pedigree = pedtools::nuclearPed(nch = 4),
ibd_state = 2, chromosome_length = 100, fraction = TRUE)
m4
stopifnot(all.equal(0.25^3, m4$mean))
Variance for IBD Distribution Objects
Description
Computes the variance of a total IBD or segment count distribution.
Usage
var(x, ...)
Arguments
x |
An object of class |
... |
Additional arguments (not used). |
Value
A numeric value.