twoLocusIBD | R Documentation |
Computes the 3*3 matrix of two-locus IBD coefficients of a pair of non-inbred pedigree members, for a given recombination rate.
twoLocusIBD(
x,
ids,
rho,
coefs = NULL,
detailed = FALSE,
uniMethod = 1,
verbose = FALSE
)
x |
A pedigree in the form of a |
ids |
A character (or coercible to character) containing ID labels of two pedigree members. |
rho |
A number in the interval |
coefs |
A character indicating which coefficient(s) to compute. A subset
of |
detailed |
A logical, indicating whether the condensed (default) or detailed coefficients should be returned. |
uniMethod |
Either 1 or 2 (for testing purposes) |
verbose |
A logical. |
Let A, B be two pedigree members, and L1, L2 two loci with a given
recombination rate \rho
. The two-locus IBD coefficients
\kappa_{i,j}(\rho)
, for 0 \le i,j \le 2
are
defined as the probability that A and B have i
alleles IBD at L1 and j
alleles IBD at L2 simultaneously. Note that IBD alleles at the two loci are
not required to be in cis (or in trans for that matter).
The method of computation depends on the (single-locus) IBD coefficient
\kappa_2
. If this is zero (e.g. if A is a direct ancestor of B, or vice
versa) the two-locus IBD coefficients are easily computable from the
two-locus kinship coefficients, as implemented in twoLocusKinship()
. In the
general case, the computation is more involved, requiring generalised
two-locus kinship coefficients. This is implemented in the function
twoLocusGeneralisedKinship()
, which is not exported yet.
By default, a symmetric 3*3 matrix containing the two-locus IBD
coefficients \kappa_{i,j}
.
If either coefs
is explicitly given (i.e., not NULL), or detailed = TRUE
, the computed coefficients are returned as a named vector.
twoLocusKinship()
, twoLocusIdentity()
, twoLocusPlot()
# Plot title used in several examples below
main = expression(paste("Two-locus IBD: ", kappa[`1,1`]))
###################################################################
# Example 1: A classic example of three relationships with the same
# one-locus IBD coefficients, but different two-locus coefficients.
# As a consequence, these relationships cannot be separated using
# unlinked markers, but are (theoretically) separable with linked
# markers.
###################################################################
peds = list(
GrandParent = list(ped = linearPed(2), ids = c(1, 5)),
HalfSib = list(ped = halfSibPed(), ids = c(4, 5)),
Uncle = list(ped = avuncularPed(), ids = c(3, 6)))
twoLocusPlot(peds, coeff = "k11", main = main, lty = 1:3, col = 1)
############################################################
# Example 2: Inspired by Fig. 3 in Thompson (1988),
# and its erratum: https://doi.org/10.1093/imammb/6.1.1.
#
# These relationships are also analysed in ?twoLocusKinship,
# where we show that they have identical two-locus kinship
# coefficients. Here we demonstrate that they have different
# two-locus IBD coefficients.
############################################################
peds = list(
GreatGrand = list(ped = linearPed(3), ids = c(1, 7)),
HalfUncle = list(ped = avuncularPed(half = TRUE), ids = c(4, 7))
)
twoLocusPlot(peds, coeff = "k11", main = main, lty = 1:2, col = 1)
########################################################################
# Example 3: Fig. 15 of Vigeland (2021).
# Two-locus IBD of two half sisters whose mother have inbreeding
# coefficient 1/4. We compare two different realisations of this:
# PO: the mother is the child of parent-offspring
# SIB: the mother is the child of full siblings
#
# The fact that these relationships have different two-locus coefficients
# exemplifies that a single-locus inbreeding coefficient cannot replace
# the genealogy in analyses of linked loci.
########################################################################
po = addChildren(nuclearPed(1, sex = 2), 1, 3, nch = 1, sex = 2)
po = addDaughter(addDaughter(po, 4), 4)
sib = addChildren(nuclearPed(2, sex = 1:2), 3, 4, nch = 1)
sib = addDaughter(addDaughter(sib, 5), 5)
# plotPedList(list(po, sib), new = TRUE, title = c("PO", "SIB"))
# List of pedigrees and ID pairs
peds = list(PO = list(ped = po, ids = leaves(po)),
SIB = list(ped = sib, ids = leaves(sib)))
twoLocusPlot(peds, coeff = "k11", main = main, lty = 1:2, col = 1)
### Check against exact formulas
rho = seq(0, 0.5, length = 11) # recombination values
kvals = sapply(peds, function(x)
sapply(rho, function(r) twoLocusIBD(x$ped, x$ids, r, coefs = "k11")))
k11.po = 1/8*(-4*rho^5 + 12*rho^4 - 16*rho^3 + 16*rho^2 - 9*rho + 5)
stopifnot(all.equal(kvals[, "PO"], k11.po, check.names = FALSE))
k11.s = 1/16*(8*rho^6 - 32*rho^5 + 58*rho^4 - 58*rho^3 + 43*rho^2 - 20*rho + 10)
stopifnot(all.equal(kvals[, "SIB"], k11.s, check.names = FALSE))
################################################
# Example 4:
# The complete two-locus IBD matrix of full sibs
################################################
x = nuclearPed(2)
k2_mat = twoLocusIBD(x, ids = 3:4, rho = 0.25)
k2_mat
### Compare with exact formulas
IBDSibs = function(rho) {
R = rho^2 + (1-rho)^2
nms = c("ibd0", "ibd1", "ibd2")
m = matrix(0, nrow = 3, ncol = 3, dimnames = list(nms, nms))
m[1,1] = m[3,3] = 0.25 *R^2
m[2,1] = m[1,2] = 0.5 * R * (1-R)
m[3,1] = m[1,3] = 0.25 * (1-R)^2
m[2,2] = 0.5 * (1 - 2 * R * (1-R))
m[3,2] = m[2,3] = 0.5 * R * (1-R)
m
}
stopifnot(all.equal(k2_mat, IBDSibs(0.25)))
#####################################################
# Example 5: Two-locus IBD of quad half first cousins
#
# We use this to exemplify two simple properties of
# the two-locus IBD matrix.
#####################################################
x = quadHalfFirstCousins()
ids = leaves(x)
# First compute the one-locus IBD coefficients (= c(17, 14, 1)/32)
k1 = kappaIBD(x, ids)
### Case 1: Complete linkage (`rho = 0`).
# In this case the two-locus IBD matrix has `k1` on the diagonal,
# and 0's everywhere else.
k2_mat_0 = twoLocusIBD(x, ids = ids, rho = 0)
stopifnot(all.equal(k2_mat_0, diag(k1), check.attributes = FALSE))
#' ### Case 2: Unlinked loci (`rho = 0.5`).
# In this case the two-locus IBD matrix is the outer product of
# `k1` with itself.
k2_mat_0.5 = twoLocusIBD(x, ids = ids, rho = 0.5)
stopifnot(all.equal(k2_mat_0.5, k1 %o% k1, check.attributes = FALSE))
########################################################
# Example 6: By Donnelly (1983) these relationships are
# genetically indistinguishable
########################################################
x1 = halfCousinPed(1)
x2 = halfCousinPed(0, removal = 2)
stopifnot(identical(
twoLocusIBD(x1, ids = leaves(x1), rho = 0.25),
twoLocusIBD(x2, ids = leaves(x2), rho = 0.25)))
###########################################################
# Example 7: Detailed coefficients of double first cousins.
# Compare with exact formulas by Denniston (1975).
###########################################################
## Not run:
x = doubleFirstCousins()
ids = leaves(x)
rho = 0.25
kapDetailed = twoLocusIBD(x, ids, rho, detailed = TRUE)
# Example 1 of Denniston (1975)
denn = function(rho) {
F = (1-rho)^2 * (rho^2 + (1-rho)^2)/4 + rho^2/8
U = 1 + 2*F
V = 1 - 4*F
# Note that some of Denniston's definitions differ slightly;
# some coefficients must be multiplied with 2
c(k00 = U^2/4,
k01 = U*V/8 *2,
k02 = V^2/16,
k10 = U*V/8 *2,
k11.cc = F*U/2 *2,
k11.ct = 0,
k11.tc = 0,
k11.tt = V^2/16 *2,
k12.h = F*V/4 *2,
k12.r = 0,
k20 = V^2/16,
k21.h = F*V/4 *2,
k21.r = 0,
k22.h = F^2,
k22.r = 0)
}
stopifnot(all.equal(kapDetailed, denn(rho)))
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.