calc_genoprob: Calculate conditional genotype probabilities

Description

Uses a hidden Markov model to calculate the probabilities of the true underlying genotypes given the observed multipoint marker data, with possible allowance for genotyping errors.

Usage

calc_genoprob(
  cross,
  map = NULL,
  error_prob = 0.0001,
  map_function = c("haldane", "kosambi", "c-f", "morgan"),
  lowmem = FALSE,
  quiet = TRUE,
  cores = 1
)

Value

An object of class "calc_genoprob": a list of three-dimensional arrays of probabilities, individuals x genotypes x positions. (Note that the arrangement is different from R/qtl.) Also contains four attributes:

crosstype - The cross type of the input cross.
is_x_chr - Logical vector indicating whether chromosomes are to be treated as the X chromosome or not, from input cross.
alleles - Vector of allele codes, from input cross.
alleleprobs - Logical value (FALSE) that indicates whether the probabilities are compressed to allele probabilities, as from genoprob_to_alleleprob().

Arguments

cross: Object of class "cross2". For details, see the R/qtl2 developer guide.
map: Genetic map of markers. May include pseudomarker locations (that is, locations that are not within the marker genotype data). If NULL, the genetic map in cross is used.
error_prob: Assumed genotyping error probability
map_function: Character string indicating the map function to use to convert genetic distances to recombination fractions.
lowmem: If FALSE, split individuals into groups with common sex and crossinfo and then precalculate the transition matrices for a chromosome; potentially a lot faster but using more memory.
quiet: If FALSE, print progress messages.
cores: Number of CPU cores to use, for parallel calculations. (If 0, use parallel::detectCores().) Alternatively, this can be links to a set of cluster sockets, as produced by parallel::makeCluster().

Details

Let \(O_k\) denote the observed marker genotype at position \(k\), and \(g_k\) denote the corresponding true underlying genotype.

We use the forward-backward equations to calculate \(\alpha_{kv} = \log Pr(O_1, \ldots, O_k, g_k = v)\) and \(\beta_{kv} = \log Pr(O_{k+1}, \ldots, O_n | g_k = v)\)

We then obtain \(Pr(g_k | O_1, \ldots, O_n) = \exp(\alpha_{kv} + \beta_{kv}) / s\) where \(s = \sum_v \exp(\alpha_{kv} + \beta_{kv})\)

Examples

Run this code

grav2 <- read_cross2(system.file("extdata", "grav2.zip", package="qtl2"))
gmap_w_pmar <- insert_pseudomarkers(grav2$gmap, step=1)
probs <- calc_genoprob(grav2, gmap_w_pmar, error_prob=0.002)