fwsim {fwsim} | R Documentation |
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 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 be the population size at generation
and
the population size at generation
.
Then we assume that
conditionally on
is
distributed for
(
gives expected growth and
gives expected decrease).
For each haplotype occuring
times in the
'th generation, the number
of children
is
distributed independently of other haplotypes.
It then follows that the sum of the number of haplotypes follows a
distribution (as just stated in the previous paragraph) and that
conditionally on
follows a
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 class 'fwsim'
print(x, ...)
## S3 method for class 'fwsim'
summary(object, ...)
## S3 method for class '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 |
N0 |
count of the |
mutmodel |
a |
alpha |
vector of length 1 or |
SNP |
to make alleles modulus 2 to immitate SNPs. |
save_generations |
to save intermediate populations. |
progress |
whether to print progress of the evolution. |
trace |
whether to print trace of the evolution (more verbose than |
ensure_children |
Ensures that every generation has at least one child; implemented by getting |
x |
A |
object |
A |
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 |
pop_sizes |
the population size for each generation |
expected_pop_sizes |
the expected population size for each generation |
Author(s)
Mikkel Meyer Andersen <mikl@math.aau.dk> and Poul Svante Eriksen
Examples
# 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)