menbayesgl {hwep}

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