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
)

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.

Arguments

x: 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.

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

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

if (FALSE) {
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)))
}