| Title: | Fisher-Wright Population Simulation |
|---|---|
| Description: | Simulates a population under the Fisher-Wright model (fixed or stochastic population size) with a one-step neutral mutation process (stepwise mutation model, logistic mutation model and exponential mutation model supported). The stochastic population sizes are random Poisson distributed and different kinds of population growth are supported. For the stepwise mutation model, it is possible to specify locus and direction specific mutation rate (in terms of upwards and downwards mutation rate). Intermediate generations can be saved in order to study e.g. drift. |
| Authors: | Mikkel Meyer Andersen and Poul Svante Eriksen |
| Maintainer: | Mikkel Meyer Andersen <[email protected]> |
| License: | GPL-2 |
| Version: | 0.3.4 |
| Built: | 2024-09-17 04:41:27 UTC |
| Source: | https://github.com/cran/fwsim |
Tools for analysing haplotypes, e.g. population simulation.
Mikkel Meyer Andersen
Maintainer: Mikkel Meyer Andersen <[email protected]>
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.
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, ...)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, ...)
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. |
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 |
Mikkel Meyer Andersen <[email protected]> and Poul Svante Eriksen
# 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)# 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)
Method to initialise a mutation model.
init_mutmodel(modeltype = 1, mutpars = NULL, ...) ## S3 method for class 'mutmodel' print(x, ...)init_mutmodel(modeltype = 1, mutpars = NULL, ...) ## S3 method for class 'mutmodel' print(x, ...)
modeltype |
1: SMM (traditional single-step mutation model). 2: LMM (Logistic mutation model introduced in Jochens (2011) 'Empirical Evaluation Reveals Best Fit of a Logistic Mutation Model for Human Y-Chromosomal Microsatellites'). 3: Exponential mutation model (unpublished). |
mutpars |
A matrix specifying the mutation parameters for each locus. Rows are parameters and columns are loci. If a vector, the same values are used for all loci. |
x |
A |
... |
not used. |
Mutation parameters for each locus.
Mutmodel 1 (SMM): 2 parameters per locus
P(i -> i-1) = mu\_d
P(i -> i+1) = mu\_u
P(i -> i) = 1 - P(i -> i-1) - P(i -> i+1)
= 1 - mu\_d - mu\_u
mutpars[1, locus]: mu\_d
mutpars[2, locus]: mu\_u
Mutmodel 2 (LMM): 6 parameters per locus
P(i -> i-1) = gamma\_d / (1 + exp(alpha\_d*(beta\_d - i)))
P(i -> i+1) = gamma\_u / (1 + exp(alpha\_u*(beta\_u - i)))
P(i -> i) = 1 - P(i -> i-1) - P(i -> i+1)
mutpars[1, locus]: gamma\_d
mutpars[2, locus]: alpha\_d
mutpars[3, locus]: beta\_d
mutpars[4, locus]: gamma\_u
mutpars[5, locus]: alpha\_u
mutpars[6, locus]: beta\_u
Mutmodel 3 (EMM): 4 parameters per locus
P(i -> i-1) = 1/((1 + exp(a + b*i))*(1 + exp(alpha + beta*i)))
P(i -> i+1) = exp(alpha+beta*i)/((1+exp(a+b*i))*(1+exp(alpha+beta*i)))
P(i -> i) = 1 - P(i -> i-1) - P(i -> i+1)
= exp(a + b*i)/(1 + exp(a + b*i))
mutpars[1, locus]: a
mutpars[2, locus]: b
mutpars[3, locus]: alpha
mutpars[4, locus]: beta
A mutmodel object (a list with entires modeltype and mutpars).
mutpars <- matrix(c(c(0.003, 0.001), rep(0.004, 2), rep(0.001, 2)), ncol = 3) colnames(mutpars) <- c("DYS19", "DYS389I", "DYS391") mutmodel <- init_mutmodel(modeltype = 1L, mutpars = mutpars) mutmodel 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) mutmodelmutpars <- matrix(c(c(0.003, 0.001), rep(0.004, 2), rep(0.001, 2)), ncol = 3) colnames(mutpars) <- c("DYS19", "DYS389I", "DYS391") mutmodel <- init_mutmodel(modeltype = 1L, mutpars = mutpars) mutmodel 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
Functions for mutation model logic, e.g. probability of downwards and upwards mutations etc.
mutmodel_not_mut(mutmodel, locus, alleles) mutmodel_dw_mut(mutmodel, locus, alleles) mutmodel_up_mut(mutmodel, locus, alleles) approx_stationary_dist(mutmodel, alleles)mutmodel_not_mut(mutmodel, locus, alleles) mutmodel_dw_mut(mutmodel, locus, alleles) mutmodel_up_mut(mutmodel, locus, alleles) approx_stationary_dist(mutmodel, alleles)
mutmodel |
a |
locus |
the locus of interest (integer, remember postfix L). |
alleles |
vector of integers (remember postfix L) of the alleles of interest. |
Mutation probabilities for locus locus in mutation model mutmodel at alleleles alleles.
Mikkel Meyer Andersen <[email protected]> and Poul Svante Eriksen
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_not_mut(mutmodel, locus = 1L, alleles = 10L:20L) mutmodel_dw_mut(mutmodel, locus = 1L, alleles = 10L:20L) mutmodel_up_mut(mutmodel, locus = 1L, alleles = 10L:20L) statdists <- approx_stationary_dist(mutmodel, alleles = 5L:20L) bp <- barplot(statdists, beside = TRUE) text(bp, 0.02, round(statdists, 1), cex = 1, pos = 3) text(bp, 0, rep(rownames(statdists), ncol(mutmodel$mutpars)), cex = 1, pos = 3) mutpars <- matrix(c(c(0.003, 0.001), rep(0.004, 2), rep(0.001, 2)), ncol = 3) colnames(mutpars) <- c("DYS19", "DYS389I", "DYS391") mutmodel <- init_mutmodel(modeltype = 1L, mutpars = mutpars) mutmodel statdists <- approx_stationary_dist(mutmodel, alleles = 5L:20L) statdists bp <- barplot(statdists, beside = TRUE) text(bp, 0.02, round(statdists, 1), cex = 1, pos = 3) text(bp, 0, rep(rownames(statdists), ncol(mutmodel$mutpars)), cex = 1, pos = 3)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_not_mut(mutmodel, locus = 1L, alleles = 10L:20L) mutmodel_dw_mut(mutmodel, locus = 1L, alleles = 10L:20L) mutmodel_up_mut(mutmodel, locus = 1L, alleles = 10L:20L) statdists <- approx_stationary_dist(mutmodel, alleles = 5L:20L) bp <- barplot(statdists, beside = TRUE) text(bp, 0.02, round(statdists, 1), cex = 1, pos = 3) text(bp, 0, rep(rownames(statdists), ncol(mutmodel$mutpars)), cex = 1, pos = 3) mutpars <- matrix(c(c(0.003, 0.001), rep(0.004, 2), rep(0.001, 2)), ncol = 3) colnames(mutpars) <- c("DYS19", "DYS389I", "DYS391") mutmodel <- init_mutmodel(modeltype = 1L, mutpars = mutpars) mutmodel statdists <- approx_stationary_dist(mutmodel, alleles = 5L:20L) statdists bp <- barplot(statdists, beside = TRUE) text(bp, 0.02, round(statdists, 1), cex = 1, pos = 3) text(bp, 0, rep(rownames(statdists), ncol(mutmodel$mutpars)), cex = 1, pos = 3)