This is the main function of the package, simulating the recombination process in each meiosis of a pedigree. The output summarises the IBD segments between all or a subset of individuals.
ibdsim(
x,
N = 1,
ids = NULL,
map = "decode",
model = c("chi", "haldane"),
skipRecomb = NULL,
simplify1 = TRUE,
seed = NULL,
verbose = TRUE
)If N > 1, 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 and in pipe chains. To enforce a list output in
this case, add simplify1 = FALSE.
A genomeSim object is essentially a numerical matrix describing the
allele flow through the pedigree in a single genome simulation. Each row
corresponds to a chromosomal segment. The first 3 columns (chrom, startMB,
endMB) describe the physical megabase location of the segment. Next come
the centimorgan 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, i.e., the number of alleles shared
identical by descent. This is always either 0, 1, 2 or NA, where the latter
appear if either individual (or both) of individuals is autozygous in the
segment.
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.
A pedtools::ped() object.
A positive integer indicating the number of simulations.
A vector of names indicating which pedigree members should be
included in the output. Alternatively, a function taking x as input and
returning such a vector, like pedtools::leaves(). By default (ids = NULL) everyone is included.
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".
Either "chi" or "haldane", indicating the statistical model for recombination (see Details). Default: "chi".
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.) Alternatively, a function taking x as
input, e.g., pedtools::founders(). By default, recombination is skipped
in founders who are uninformative for IBD sharing in the ids individuals.
A logical, by default TRUE, removing the outer list layer when N = 1. See Value.
An integer to be passed on to set.seed()).
A logical.
Each simulation starts by distributing unique alleles (labelled 1, 2, ...) to
the pedigree founders, followed by a separate recombination process for each
chromosome, according to the value of model:
model = "chi" (default): This uses a renewal process along the
four-strand bundle, with waiting times following a chi-square distribution.
model = "haldane": In this model, crossover events are modelled as a
Poisson process along each chromosome. (Faster, but less realistic.)
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 may be
replaced by simpler maps with constant recombination rates. See
uniformMap() and loadMap() for ways to produce such maps.
Halldorsson et al. Characterizing mutagenic effects of recombination through a sequence-level genetic map. Science 363, no. 6425 (2019).
### 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))
Run the code above in your browser using DataLab