Simulator details

In this package, synthetic sequencing data are generated using a statistical model for an individual in a given population. Each individual is assigned a COI value, which is used to simulate the number of sequence reads mapped to the reference and alternative allele for the biallelic SNPs considered. These are then used to derive the within-sample frequency of the population-level minor allele at each locus i, denoted as w_i.

We first define the number of loci being simulated, l. We assume that the distribution of reference allele frequencies for locus i, R_i, is described by a Beta distribution with shape parameters α and β: $$ \mathbb{P}(R_i = r_i) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}r_i^{(\alpha-1)}(1-r_i)^{(\beta-1)}, $$

where Γ is the Gamma function. In our simulations, we assume α = 1 and β = 5 and draw l values of R_i. In our methods, we have defined relationships with respect to the frequency of the minor allele, p_i. Recall that p_i ∈ [0,0.5]. Thus, we define $$ p_i = \begin{cases} r_i & r_i \leq 0.5 \\ 1 - r_i & r_i > 0.5. \end{cases} $$

To simulate w_i for each individual, we first assign the COI, k. Next, we simulate the genotype at locus i, drawing k values from a Binomial distribution with probability p_i, i.e., the probability of having the minor allele at locus i is equal to the population-level frequency of the minor allele. This yields the matrix, 𝔾 ∈ ℝ^k × l, with k rows and l columns defining the phased haplotype of each strain, where 𝔾_i, j = 1 indicates that strain j at locus i has the minor allele.

To determine the number of sequence reads related to each strain, we first draw the proportion of each strain, S = (S₁,…,S_k), in an individual with k strains, which we assume is described by a Dirichlet distribution with concentration parameters α₁, …, α_k. The “true” WSMAF of the minor allele at locus i, denoted τ_i, is thus given by the sum of the strain proportions for strains with the minor allele. We can express this using matrix multiplication: τ = S × 𝔾_i, where τ ∈ ℝ^1 × l, S ∈ ℝ^1 × k, and 𝔾 ∈ ℝ^k × l.

In simulations with no assumed sequence error, the number of sequence reads with the minor allele at each locus is given by sampling from a Binomial distribution c times with probability τ, where c is the assumed number of sequence reads, i.e., the read depth or coverage. However, in simulations with sequence error, we perturb the “true” WSMAF by assuming that a number of sequence reads are incorrectly sampled. An incorrectly sampled locus with the major allele will yield the minor allele and vice versa. In our simulations, we assume a fixed sequence error, ϵ, such that the probability of correctly sampling the minor allele is 1 − ϵ, and the probability of incorrectly sampling the minor allele is ϵ. Therefore, we can represent the probability of sampling the minor allele, now denoted as τ_ϵ to indicate the added error, as τ_ϵ = τ(1−ϵ) + (1−τ)ϵ.

Lastly, we may account for overdispersion—additional unexpected variability in our data—and sample the number of sequence reads with the minor allele at each locus from a Beta-binomial distribution rather than a Binomial distribution. Dividing the number of instances of the minor allele by the coverage at each locus results in the simulated WSMAF.

simulation <- sim_biallelic(2, plmaf = runif(1000, 0, 0.5))
plot(simulation)

# Set the seed
set.seed(1)

# Define number of loci and the PLMAF
n_loci <- 5
plmaf <- rbeta(n_loci, 1, 5)
plmaf[plmaf > 0.5] <- 1 - plmaf[plmaf > 0.5]

strain_p <- coiaf:::rdirichlet(rep(alpha, COI))

haplotype <- mapply(function(x) rbinom(COI, 1, x), x = plmaf)

true_wsmaf <- colSums(sweep(haplotype, 1, strain_p, "*"))

# Rounding errors can cause numbers greater than 1
true_wsmaf[true_wsmaf > 1] <- 1L

# Account for sequencing error
true_wsmaf <- true_wsmaf * (1 - epsilon) + (1 - true_wsmaf) * epsilon

# Sample sequence reads
counts <- rbinom(n_loci, size = rep(coverage, n_loci), prob = true_wsmaf)

# Determine simulated WSMAF
sim_wsmaf <- counts / coverage