ibdsim {ibdsim2}

IBD simulation

Description

This is the main function of the package, simulating the recombination process in each meioses of a pedigree. The output summarises the IBD segments between all or a subset of individuals.

Usage

ibdsim(
  x,
  N = 1,
  ids = labels(x),
  map = "decode",
  model = c("chi", "haldane"),
  skipRecomb = NULL,
  simplify1 = TRUE,
  seed = NULL,
  verbose = TRUE
)

Arguments

x

A pedtools::ped() object.

N

A positive integer indicating the number of simulations.

ids

A subset of pedigree members whose IBD sharing should be analysed. If NULL, all members are included.

map

The genetic map to be used in the simulations: Allowed values are:

• a genomeMap object, typically produced by loadMap()

• a single chromMap object, for instance as produced by uniformMap()

• a character, which is passed on to loadMap() with default parameters. Currently the only valid option is "decode19" (or abbreviations of this).

Default: "decode19".

model

Either "chi" or "haldane", indicating the statistical model for recombination (see details). Default: "chi".

skipRecomb

A vector of ID labels indicating individuals whose meioses should be simulated without recombination. (Each child will then receive a random strand of each chromosome.) The default action is to skip recombination in founders who are uninformative for IBD sharing in the ids individuals.

simplify1

A logical, by default TRUE, removing the outer list layer when N = 1. See Value.

seed

An integer to be passed on to set.seed()).

verbose

A logical.

Details

Each simulation starts by unique alleles (labelled 1, 2, ...) being distributed to the pedigree founders. In each meiosis, homologue chromosomes are made to recombine according to the value of model:

• model = "haldane": In this model, crossover events are modelled as a Poisson process along each chromosome.

• model = "chi" (default): This uses a renewal process along the four-strand bundle, with waiting times following a chi square distribution.

Recombination rates along each chromosome are determined by the map parameter. The default value ("decode19") loads a thinned version of the recombination map of the human genome published by Halldorsson et al (2019).

In many applications, the fine-scale default map is not necessary, and should be replaced by simpler maps with constant recombination rates. See uniformMap() and loadMap() for ways to produce such maps.

Value

A list of N objects of class genomeSim.

If N = 1 the outer list layer is removed by default, which is typically desired in interactive use (especially when piping). To enforce a list output, add simplify1 = FALSE.

A genomeSim object is essentially a numerical matrix describing the allele flow through the pedigree in a single simulated. Each row corresponds to a chromosomal segment. The first 3 columns (chrom, startMB, endMB) describe the physical location of the segment. Next, the genetic coordinates (startCM, endCM), which are computed from map by averaging the male and female values. Then follow the allele columns, two for each individual in ids, suffixed by ":p" and ":m" signifying the paternal and maternal alleles, respectively.

If ids has length 1, a column named Aut is added, whose entries are 1 for autozygous segments and 0 otherwise.

If ids has length 2, two columns are added:

• IBD : The IBD status of each segment (= number of alleles shared identical by descent). For a given segment, the IBD status is either 0, 1, 2 or NA. If either individual is autozygous in a segment, the IBD status is reported as NA. With inbred individuals the Sigma column (see below) is more informative than the IBD column.

• Sigma : The condensed identity ("Jacquard") state of each segment, given as an integer in the range 1-9. The numbers correspond to the standard ordering of the condensed states. In particular, for non-inbred individuals, the states 9, 8, 7 correspond to IBD status 0, 1, 2 respectively.

References

Halldorsson et al. Characterizing mutagenic effects of recombination through a sequence-level genetic map. Science 363, no. 6425 (2019).

Examples

### Example 1: Half siblings ###

hs = halfSibPed()
sim = ibdsim(hs, map = uniformMap(M = 1), seed = 10)
sim

# Plot haplotypes
haploDraw(hs, sim)

#' ### Example 2: Full sib mating ###

x = fullSibMating(1)
sim = ibdsim(x, ids = 5:6, map = uniformMap(M = 10), seed = 1)
head(sim)

# All 9 identity states are present
stopifnot(setequal(sim[, 'Sigma'], 1:9))