menbayesgl: Bayes test for F1/S1 genotype frequencies using genotype...

In hwep: Hardy-Weinberg Equilibrium in Polyploids

Bayes test for F1/S1 genotype frequencies using genotype likelihoods

Description

Uses get_q_array() from the updog R package to calculate segregation probabilities (assuming no double reduction) and tests that offspring genotypes follow this distribution.

Usage

menbayesgl(
  gl,
  method = c("f1", "s1"),
  p1gl = NULL,
  p2gl = NULL,
  lg = TRUE,
  beta = NULL,
  chains = 2,
  cores = 1,
  iter = 2000,
  warmup = floor(iter/2),
  ...
)

Arguments

gl	A matrix of genotype log-likelihoods. The rows index the individuals and the columns index the genotypes. So gl[i,k] is the genotype log-likelihood for individual i at dosage k-1. We assume the natural log is used (base e).
method	Should we test for F1 proportions ("f1") or S1 proportions ("s1")?
p1gl	A vector of genotype log-likelihoods for parent 1. p1gl[k] is the log-likelihood of parent 1's data given their genotype is k.
p2gl	A vector of genotype log-likelihoods for parent 2. p2gl[k] is the log-likelihood of parent 2's data given their genotype is k.
lg	A logical. Should we return the log Bayes factor (TRUE) or the Bayes factor (FALSE)?
beta	The concentration hyperparameters of the genotype frequencies under the alternative of no random mating. Should be length ploidy + 1.
chains	Number of MCMC chains. Almost always 1 is enough, but we use 2 as a default to be conservative.
cores	Number of cores to use.
iter	Total number of iterations.
warmup	Number of those iterations used in the burnin.
...	Control arguments passed to sampling().

Author(s)

David Gerard

References

• Gerard D (2022). "Bayesian tests for random mating in autopolyploids." bioRxiv. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1101/2022.08.11.503635")}.

Examples

## Not run: 
set.seed(1)
ploidy <- 4

## Simulate under the null ----
q <- updog::get_q_array(ploidy = 4)[3, 3, ]

## See BF increases
nvec <- c(stats::rmultinom(n = 1, size = 10, prob = q))
gl <- simgl(nvec = nvec)
menbayesgl(gl = gl, method = "f1")

nvec <- c(stats::rmultinom(n = 1, size = 100, prob = q))
gl <- simgl(nvec = nvec)
menbayesgl(gl = gl, method = "f1")

nvec <- c(stats::rmultinom(n = 1, size = 1000, prob = q))
gl <- simgl(nvec = nvec)
menbayesgl(gl = gl, method = "f1")

## Simulate under the alternative ----
q <- stats::runif(ploidy + 1)
q <- q / sum(q)

## See BF decreases
nvec <- c(stats::rmultinom(n = 1, size = 10, prob = q))
gl <- simgl(nvec = nvec)
menbayesgl(gl = gl, method = "f1")

nvec <- c(stats::rmultinom(n = 1, size = 100, prob = q))
gl <- simgl(nvec = nvec)
menbayesgl(gl = gl, method = "f1")

nvec <- c(stats::rmultinom(n = 1, size = 1000, prob = q))
gl <- simgl(nvec = nvec)
menbayesgl(gl = gl, method = "f1")


## End(Not run)

hwep documentation built on May 31, 2023, 9:06 p.m.