fwsim (version 0.3.4)

fwsim: Fisher-Wright Population Simulation

Description

This package provides tools to simulate a population under the Fisher-Wright model with a stepwise neutral mutation process on \(r\) loci, where mutations on loci happen independently. The population sizes are either fixed (traditional/original Fisher-Wright model) or random Poisson distributed with exponential growth supported. Intermediate generations can be saved in order to study e.g. drift.

For stochastic population sizes: Model described in detail at http://arxiv.org/abs/1210.1773. Let \(M\) be the population size at generation \(i\) and \(N\) the population size at generation \(i + 1\). Then we assume that \(N\) conditionally on \(M\) is \(\mbox{Poisson}(\alpha M)\) distributed for \(\alpha > 0\) (\(\alpha > 1\) gives expected growth and \(0 < \alpha < 1\) gives expected decrease).

For each haplotype \(x\) occuring \(m\) times in the \(i\)'th generation, the number of children \(n\) is \(\mbox{Poisson}(\alpha m)\) distributed independently of other haplotypes. It then follows that the sum of the number of haplotypes follows a \(\mbox{Poisson}(\alpha M)\) distribution (as just stated in the previous paragraph) and that \(n\) conditionally on \(N\) follows a \(\mbox{Binomial}(N, m/M)\) as expected.

The mutation model can be e.g. the stepwise neutral mutation model. See init_mutmodel for details.

Usage

fwsim(G, H0, N0, mutmodel, alpha = 1.0, SNP = FALSE, 
  save_generations = NULL, progress = TRUE, trace = FALSE, ensure_children = FALSE, ...)
fwsim_fixed(G, H0, N0, mutmodel, SNP = FALSE, 
  save_generations = NULL, progress = TRUE, trace = FALSE, ...)
# S3 method for fwsim
print(x, ...)
# S3 method for fwsim
summary(object, ...)
# S3 method for fwsim
plot(x, which = 1L, ...)

Arguments

G

number of generations to evolve (integer, remember postfix L).

H0

haplotypes of the initial population. Must be a vector or matrix (if more than one initial haplotype). The number of loci is the length or number of columns of H0.

N0

count of the H0 haplotypes. The i'th element is the count of the haplotype H0[i, ]. sum(N0) is the size of initial population.

mutmodel

a mutmodel object created with init_mutmodel. Alternatively, a numeric vector of length r of mutation probabilities (this will create a stepwise mutation model with r loci divide the mutation probabilities evenly between upwards and downwards mutation).

alpha

vector of length 1 or G of growth factors (1 correspond to expected constant population size). If length 1, the value is reused in creating a vector of length G.

SNP

to make alleles modulus 2 to immitate SNPs.

save_generations

to save intermediate populations. NULL means that no intermediate population will be saved. Else, a vector of the generation numbers to save.

progress

whether to print progress of the evolution.

trace

whether to print trace of the evolution (more verbose than progress).

ensure_children

Ensures that every generation has at least one child; implemented by getting \(\mbox{Poisson}(\alpha M) + 1\) children.

x

A fwsim object.

object

A fwsim object.

which

A number specifying the plot (currently only 1: the actual population sizes vs the expected sizes).

…

not used.

Value

A fwsim object with elements

pars

the parameters used for the simulation

saved_populations

a list of haplotypes in the intermediate populations

population

haplotypes in the end population after G generations

pop_sizes

the population size for each generation

expected_pop_sizes

the expected population size for each generation

Examples

# NOT RUN {
# SMM (stepwise mutation model) example
set.seed(1)
fit <- fwsim(G = 100L, H0 = c(0L, 0L, 0L), N0 = 10000L, 
  mutmodel = c(Loc1 = 0.001, Loc2 = 0.002, Loc3 = 0.003))
summary(fit)
fit

# SMM (stepwise mutation model) example
H0 <- matrix(c(0L, 0L, 0L), 1L, 3L, byrow = TRUE)

mutmodel <- init_mutmodel(modeltype = 1L, 
  mutpars = matrix(c(c(0.003, 0.001), rep(0.004, 2), rep(0.001, 2)), 
                   ncol = 3, 
                   dimnames = list(NULL, c("DYS19", "DYS389I", "DYS391"))))
mutmodel

set.seed(1)
fit <- fwsim(G = 100L, H0 = H0, N0 = 10000L, mutmodel = mutmodel)

xtabs(N ~ DYS19 + DYS389I, fit$population)
plot(1L:fit$pars$G, fit$pop_sizes, type = "l", 
  ylim = range(range(fit$pop_sizes), range(fit$expected_pop_sizes)))
points(1L:fit$pars$G, fit$expected_pop_sizes, type = "l", col = "red")

set.seed(1)
fit_fixed <- fwsim_fixed(G = 100L, H0 = H0, N0 = 10000L, mutmodel = mutmodel)

# LMM (logistic mutation model) example
mutpars.locus1 <- c(0.149,   2.08,    18.3,   0.149,   0.374,   27.4) # DYS19
mutpars.locus2 <- c(0.500,   1.18,    18.0,   0.500,   0.0183,  349)  # DYS389I
mutpars.locus3 <- c(0.0163,  17.7,    11.1,   0.0163,  0.592,   14.1) # DYS391
mutpars <- matrix(c(mutpars.locus1, mutpars.locus2, mutpars.locus3), ncol = 3)
colnames(mutpars) <- c("DYS19", "DYS389I", "DYS391")
mutmodel <- init_mutmodel(modeltype = 2L, mutpars = mutpars)
mutmodel

set.seed(1)
H0_LMM <- matrix(c(15L, 13L, 10L), 1L, 3L, byrow = TRUE)
fit_LMM <- fwsim(G = 100L, H0 = H0_LMM, N0 = 10000L, mutmodel = mutmodel)
xtabs(N ~ DYS19 + DYS389I, fit_LMM$population)
# }