Bayes test for F1/S1 genotype frequencies using genotype likelihoods

Source: R/bayes.R

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.

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

References

  • Gerard D (2023). "Bayesian tests for random mating in polyploids." Molecular Ecology Resources, In press. doi:10.1111/1755-0998.13856 .

Author

David Gerard

Examples

if (FALSE) {
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")

}