Simulate biallelic data

Simulate biallelic data from a simple statistical model. Inputs include the complexity of infection and some parameters dictating skew and error distributions. Outputs include the phased haplotypes and the un-phased read count and coverage data.

sim_biallelic(
  COIs = c(1, 1),
  haplotypematrix = matrix(1, 2, 2),
  shape1 = 1,
  shape2 = 1,
  coverage = 100,
  alpha = 1,
  overdispersion = 0,
  epsilon = 0
)

Arguments

COIs

integer vector of the number of haplotypes within each sample.

haplotypematrix

matrix of haplotypes that correspond to within individual COI.

shape1

the alpha value for the beta binomial distribution, which the population-allele frequency for each allele is drawn from.

shape2

the beta value for the beta binomial distribution, which the population-allele frequency for each allele is drawn from.

coverage

coverage at each locus. If a single value then the same coverage is applied over all loci.

alpha

shape parameter of the symmetric Dirichlet prior on strain proportions.

overdispersion

the extent to which counts are over-dispersed relative to the binomial distribution. Counts are Beta-binomially distributed, with the beta distribution having shape parameters p/overdispersion and (1-p)/overdispersion.

epsilon

the probability of a single read being mis-called as the other allele. Applies in both directions.

Value

List of non-referent within-sample allele frequency dataframe (unphased) and phased vectors of non-referent within-sample allele counts, overall coverage, strain proportions, and the biallelic haplotype matrix.

Details

Simulated data are drawn from a simple statistical model:

1. Strain proportions are drawn from a symmetric Dirichlet distribution with shape parameter alpha.

2. Phased haplotypes are drawn at every locus, one for each COI. The allele at each locus is drawn from a Bernoulli distribution with probability given by the PLAF.

3. The "true" within-sample allele frequency at every locus is obtained by multiplying haplotypes by their strain proportions, and summing over haplotypes. Errors are introduced through the equation $$wsaf_{error} = wsaf*(1-e) + (1-wsaf)*e$$ where \(wsaf\) is the WSAF without error and \(e\) is the error parameter epsilon.

4. Final read counts are drawn from a beta-binomial distribution with expectation \(w_{error}\). The raw number of draws is given by the coverage, and the skew of the distribution is given by the overdispersion parameter. If overdispersion = 0 then the distribution is binomial, rather than beta-binomial.