SFS simulation

SFS simulation#

class Simulation(sfs_neut: Spectrum, params: dict = None, model: Parametrization | str = 'GammaExpParametrization', intervals_del: Tuple[float, float, int] = (-100000000.0, -1e-05, 1000), intervals_ben: Tuple[float, float, int] = (1e-05, 10000.0, 1000), h_callback: Callable[[ndarray], ndarray] = None, integration_mode: Literal['midpoint', 'quad'] = 'midpoint', linearized: bool = True, discretization: Discretization = None, parallelize: bool = True)[source]#

Bases: object

Simulate an SFS under selection given a DFE and a neutral SFS.

Example usage:

import fastdfe as fd

# create simulation object by specifying neutral SFS and DFE
sim = fd.Simulation(
    sfs_neut=fd.Simulation.get_neutral_sfs(n=20, n_sites=1e8, theta=1e-4),
    params=dict(S_d=-300, b=0.3, p_b=0.1, S_b=0.1),
    model=fd.GammaExpParametrization()
)

# perform the simulation
sim.run()

# plot spectra
sim.get_spectra().plot()
__init__(sfs_neut: Spectrum, params: dict = None, model: Parametrization | str = 'GammaExpParametrization', intervals_del: Tuple[float, float, int] = (-100000000.0, -1e-05, 1000), intervals_ben: Tuple[float, float, int] = (1e-05, 10000.0, 1000), h_callback: Callable[[ndarray], ndarray] = None, integration_mode: Literal['midpoint', 'quad'] = 'midpoint', linearized: bool = True, discretization: Discretization = None, parallelize: bool = True)[source]#

Create a simulation object.

Parameters:

  • sfs_neut (Spectrum) – Neutral SFS. This sfs is informative on the population sample size, population mutation rate, the number of sites, and demography. get_neutral_sfs() can be used to obtain a neutral SFS.

  • params (dict) – Parameters of the DFE parametrization (see model) in addition to ancestral misidentification params eps and dominance h. By default, eps=0, h=0.5, the default parameters of model will be used.

  • model (Parametrization | str) – DFE parametrization model

  • intervals_del (Tuple[float, float, int]) – (start, stop, n_interval) for deleterious population-scaled selection coefficients. The intervals will be log10-spaced. Decreasing the number of intervals to 100 provides nearly identical results while increasing speed, especially when precomputing across dominance coefficients.

  • intervals_ben (Tuple[float, float, int]) – Same as intervals_del but for positive selection coefficients. Decreasing the number of intervals to 100 provides nearly identical results while increasing speed, especially when precomputing across dominance coefficients.

  • h_callback (Callable[[ndarray], ndarray]) – A function mapping the scalar parameter h and the array of selection coefficients S to dominance coefficients of the same shape, allowing models where h depends on S. The default is lambda h, S: np.full_like(S, h), keeping h constant. Expected allele counts for a given dominance value are obtained by linear interpolation between precomputed values in intervals_h. The inferred parameter is still named h, even if transformed by h_callback, and its bounds, scales, and initial values can be set via bounds, scales, and x0. The fitness of heterozygotes and mutation homozygotes is defined as 1 + 2hs and 1 + 2s, respectively.

  • integration_mode (Literal['midpoint', 'quad']) – Integration mode when computing expected SFS under semidominance. quad is not recommended.

  • linearized (bool) – Whether to discretize and cache the linearized integral mapping DFE to SFS or use scipy.integrate.quad in each call. False not recommended.

  • parallelize (bool) – Whether to parallelize computations

model: Parametrization#

The DFE parametrization

params: dict#

Parameters for the DFE parametrization

sfs_neut: Spectrum#

Neutral SFS

sfs_sel: Spectrum#

Modelled SFS

dfe: DFE#

DFE instance

n: int#

SFS sample size

theta: float#

Population mutation rate

n_sites: float#

Number of sites

discretization: Discretization#

Discretization instance

run()[source]#

Simulate an SFS given a DFE.

Return type:

Spectrum

Returns:

Simulated SFS

get_spectra()[source]#

Get the spectra used in the simulation.

Return type:

Spectra

Returns:

Tuple of neutral SFS and selected SFS

get_alpha()[source]#

Get the proportion of adaptive substitutions (alpha) under the specified DFE and parameters.

Return type:

float

Returns:

Proportion of adaptive substitutions

static get_neutral_sfs(theta: float, n_sites: float, n: int, r: Sequence[float] = None)[source]#

Obtain a standard neutral SFS for a given theta and number of sites.

Parameters:

  • theta (float) – Population mutation rate

  • n_sites (float) – Number of total sites

  • n (int) – Number of frequency classes

  • r (Sequence[float]) – Nuisance parameters that account for demography. An array of length n-1 whose elements are multiplied element-wise with the polymorphic counts of the Kingman SFS. By default, no demography effects are considered which is equivalent to r = [1] * (n-1). Note that non-default values of r will also affect estimates of the population mutation rate.

Return type:

Spectrum

Returns:

Neutral SFS

get_wright_fisher(pop_size: int, generations: int)[source]#

Get a Wright-Fisher simulation object.

Parameters:

  • pop_size (int) – Effective population size

  • generations (int) – Number of generations to simulate

Return type:

WrightFisherSimulation

Returns:

Wright-Fisher simulation object

Contents