twoLocusIBD: Two-locus IBD coefficients

Description

Computes the 3*3 matrix of two-locus IBD coefficients of a pair of non-inbred pedigree members, for a given recombination rate.

Usage

twoLocusIBD(
  x,
  ids,
  rho,
  coefs = NULL,
  detailed = FALSE,
  uniMethod = 1,
  verbose = FALSE
)

Arguments

A pedigree in the form of a pedtools::ped object.

ids

A character (or coercible to character) containing ID labels of two pedigree members.

rho

A number in the interval \([0, 0.5]\); the recombination rate between the two loci.

coefs

A character indicating which coefficient(s) to compute. A subset of c('k00', 'k01', 'k02', 'k10', 'k11', 'k12', 'k20', 'k21', 'k22'). By default, all coefficients are computed.

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.

Value

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.

Details

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.

Examples

Run this code

# NOT RUN {
# Some variables used in several examples below
rseq = seq(0, 0.5, length = 11)  # recombination values

xlab = "Recombination rate"
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 = cousinPed(0, 1), ids = c(3, 6)))

# Compute `k11` for each rho
kvals = sapply(peds, function(x)
  sapply(rseq, function(r) twoLocusIBD(x$ped, x$ids, r, coefs = "k11")))

# Plot
matplot(rseq, kvals, type = "l", xlab = xlab, ylab = "", main = main)
legend("topright", names(peds), col = 1:3, lty = 1:3)


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

# List of pedigrees and ID pairs
peds = list(
  GreatGrand = list(ped = linearPed(3),        ids = c(1, 7)),
  HalfUncle  = list(ped = halfCousinPed(0, 1), ids = c(3, 7))
)

# Compute `k11` for each rho
kvals = sapply(peds, function(x)
  sapply(rseq, function(r) twoLocusIBD(x$ped, x$ids, r, coefs = "k11")))

# Plot
matplot(rseq, kvals, type = "l", xlab = xlab, ylab = "", main = main)
legend("topright", names(peds), col = 1:2, lty = 1:2)


######################################################################
# Example 3: 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
#
# We show below that these relationships have different two-locus
# coefficients. This exemplifies that a single-locus inbreeding
# coefficient cannot replace the genealogy in analyses of linked loci.
######################################################################

xPO = addChildren(nuclearPed(1, sex = 2), 1, 3, nch = 1, sex = 2)
xPO = addDaughter(addDaughter(xPO, 4), 4)

xSIB = addChildren(nuclearPed(2, sex = 1:2), 3, 4, nch = 1)
xSIB = addDaughter(addDaughter(xSIB, 5), 5)

plotPedList(list(xPO, xSIB), new = TRUE, title = c("PO", "SIB"))

# List of pedigrees and ID pairs
peds = list(PO = list(ped = xPO, ids = c(6, 8)),
            SIB = list(ped = xSIB, ids = c(7, 9)))

# Compute `k11` for each rho
kvals = sapply(peds, function(x)
  sapply(rseq, function(r) twoLocusIBD(x$ped, x$ids, r, coefs = "k11")))

# Plot
matplot(rseq, kvals, type = "l", xlab = xlab, ylab = "", main = main)
legend("topright", names(peds), col = 1:2, lty = 1:2)

# Check against exact formula
r = rseq
k11_PO = 1/8*(-4*r^5 + 12*r^4 - 16*r^3 + 16*r^2 - 9*r + 5)
stopifnot(all.equal(kvals[, "PO"], k11_PO, check.names = FALSE))

k11_S = 1/16*(8*r^6 - 32*r^5 + 58*r^4 - 58*r^3 + 43*r^2 - 20*r + 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 explicit 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 = c(9, 10)

# 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)))


# }