Distributions of effect sizes#
Note
The lists of objects described here are passed to the sregions parameter when initializing instances of fwdpy11.ModelParams.
One of the main reasons to perform forward simulations is to be able to model mutations affecting individual fitness. To do so, we need to specify both mutation rates and the resulting effect sizes.
fwdpy11 works by specifying an overall mutation rate to variants affecting fitness.
Given that a mutation occurs, we need to specify its “effect size”.
fwdpy11 chooses the effect size of a new mutation by first determining what region is mutated and then generating a mutation from the distribution of effect size associated with that region.
Each region is represented by instances of classes derived from the ABC fwdpy11.Sregion.
Each instance is associated with a weight.
These weights establish the relative probability that a mutation comes from a given region.
Thus, given an overall mutation rate to non-neutral variants, instances of “sregions” are used to set up a multinomial distribution for generating new mutations.
The following sets up a model where mutations have a constant effect size (\(s=-0.01\)), dominance \(h=0.25\), and occur uniformly on the interval \([0, 1)\):
import fwdpy11
sregions = [fwdpy11.ConstantS(beg=0.0, end=1.0, weight=1.0, s=-0.01, h=0.25)]
The previous example uses argument names for clarity, and the following is equivalent, with the int values getting converted to float automatically:
sregions = [fwdpy11.ConstantS(0, 1, 1, -0.01, 0.25)]
sregions[0]
fwdpy11.ConstantS(beg=0, end=1, weight=1, s=-0.01, h=0.25, coupled=True, label=0, scaling=1.0)
Note that the constructor parameters for these classes often have default values–see the specific class documentation for details.
In some scenarios, it is useful to think about the distribution of effect sizes as scaled with respect to the population size.
For example, selection coefficients may be exponentially-distributed with a mean of \(2Ns\).
To do this in fwdpy11:
# ALPHA = 2Ns
MEAN_ALPHA = -10
N = 1000
sregions = [fwdpy11.ExpS(0, 1, 1, MEAN_ALPHA, scaling=2 * N)]
sregions[0]
fwdpy11.ExpS(beg=0, end=1, weight=1, mean=-10, h=1.0, coupled=True, label=0, scaling=2000)
Region weighting#
When multiple “sregion” objects are used, the default behavior is to multiply the input weight by end-beg:
sregions = [
fwdpy11.ExpS(beg=0.0, end=1.0, weight=1.0, mean=-0.2),
fwdpy11.ConstantS(beg=1.0, end=3.0, weight=1.0, s=-0.1),
]
sregions
[fwdpy11.ExpS(beg=0.0, end=1.0, weight=1.0, mean=-0.2, h=1.0, coupled=True, label=0, scaling=1.0),
fwdpy11.ConstantS(beg=1.0, end=3.0, weight=2.0, s=-0.1, h=1.0, coupled=True, label=0, scaling=1.0)]
Here, the input weight is interpreted to mean the weight “per site” is constant.
In this example, twice as many mutations will have positions in \([1, 3)\) as from \([0, 1)\).
To change the default behavior, one can prevent the coupling between input weight and region length:
sregions = [
fwdpy11.ExpS(beg=0.0, end=1.0, weight=1.0, mean=-0.2, coupled=False),
fwdpy11.ConstantS(beg=1.0, end=3.0, weight=1.0, s=-0.1, coupled=False),
]
sregions
[fwdpy11.ExpS(beg=0.0, end=1.0, weight=1.0, mean=-0.2, h=1.0, coupled=False, label=0, scaling=1.0),
fwdpy11.ConstantS(beg=1.0, end=3.0, weight=1.0, s=-0.1, h=1.0, coupled=False, label=0, scaling=1.0)]
The absolute values of the weight parameters themselves is irrelevant.
The only thing that matters is the relative values from region to region.
Simulations based on the above examples would give the same results if the weight were 42 or 73.06.
Therefore, we can recreate our first example with code like the following:
sregions = [
fwdpy11.ExpS(beg=0.0, end=1.0, weight=56.0, mean=-0.2, coupled=False),
fwdpy11.ConstantS(beg=1.0, end=3.0, weight=112.0, s=-0.1, coupled=False),
]
sregions
[fwdpy11.ExpS(beg=0.0, end=1.0, weight=56.0, mean=-0.2, h=1.0, coupled=False, label=0, scaling=1.0),
fwdpy11.ConstantS(beg=1.0, end=3.0, weight=112.0, s=-0.1, h=1.0, coupled=False, label=0, scaling=1.0)]
In the above example, twice as many mutations occur in the second region because the weights have relative values of 2:1.
Note
Different regions are allowed to overlap, allowing the simulation of concepts like “coding regions” where the DFE are a weighted mixture from multiple distributions, etc.
Setting the dominance of new mutations#
The dominance of a new mutation is set by the h parameter during initialization:
fwdpy11.ExpS(beg=0.0, end=1.0, weight=1.0, mean=-0.2, h=1.0)
fwdpy11.ExpS(beg=0.0, end=1.0, weight=1.0, mean=-0.2, h=0.0)
fwdpy11.ExpS(beg=0.0, end=1.0, weight=1.0, mean=-0.2, h=0.0, coupled=True, label=0, scaling=1.0)
Built-in distributions of effect sizes#
Labelling mutations from different regions#
It may be of use to know what region a mutation came from.
To do so, give a nonzero value to the label argument:
fwdpy11.ConstantS(beg=0.0, end=1.0, weight=1.0, s=1e-3, label=1)
fwdpy11.ConstantS(beg=0.0, end=1.0, weight=1.0, s=0.001, h=1.0, coupled=True, label=1, scaling=1.0)
At the end of the simulation, mutations from this region will have the label value stored in the attribute fwdpy11.Mutation.label.
The value of label must fit into a 16-bit unsigned integer, e.g., numpy.uint16.
Larger values, or negative values, will result in exceptions.
The following example tries to use a value one larger than the maximum allowed:
import numpy as np
fwdpy11.ConstantS(
beg=0.0, end=1.0, weight=1.0, s=1e-3, label=np.iinfo(np.uint16).max + 1
)
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[9], line 3
1 import numpy as np
----> 3 fwdpy11.ConstantS(
4 beg=0.0, end=1.0, weight=1.0, s=1e-3, label=np.iinfo(np.uint16).max + 1
5 )
File <attrs generated init fwdpy11.regions.ConstantS>:12, in __init__(self, beg, end, weight, s, h, coupled, label, scaling)
10 _inst_dict['label'] = label
11 _inst_dict['scaling'] = scaling
---> 12 self.__attrs_post_init__()
File ~/work/fwdpy11/fwdpy11/fwdpy11/regions.py:136, in ConstantS.__attrs_post_init__(self)
135 def __attrs_post_init__(self):
--> 136 super(ConstantS, self).__init__(
137 self.beg,
138 self.end,
139 self.weight,
140 self.s,
141 self.h,
142 self.coupled,
143 self.label,
144 self.scaling,
145 )
TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
1. fwdpy11._fwdpy11._ll_ConstantS(beg: float, end: float, weight: float, s: float, h: float, coupled: bool, label: int, scaling: float)
2. fwdpy11._fwdpy11._ll_ConstantS(beg: float, end: float, weight: float, s: float, h: fwdpy11._fwdpy11.MutationDominance, coupled: bool, label: int, scaling: float)
Invoked with: 0.0, 1.0, 1.0, 0.001, 1.0, True, 65536, 1.0