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