Tracking user-specified new mutations#

Hide code cell source 
import fwdpy11
import numpy as np
import msprime

import fwdpy11.conditional_models
Hide code cell source 
def setup(prune_selected=False):
    # Dropping mutations requires existing
    # ancestry, which we can get either
    # from a burn-in or from msprime.
    initial_ts = msprime.sim_ancestry(
        samples=500,
        population_size=500,
        recombination_rate=1e-1,
        random_seed=43215,
        sequence_length=1.0,
    )

    # Build the pop from msprime output
    pop = fwdpy11.DiploidPopulation.create_from_tskit(initial_ts)

    # Set up basic model parameters
    pdict = {
        "recregions": [fwdpy11.PoissonInterval(0, 1, 1e-1)],
        # Here, 2 means that fitness is multiplicative
        # over 1, 1+hs, 1+2s.
        "gvalue": fwdpy11.Multiplicative(2.0),
        "rates": (0, 0, None),
        "prune_selected": False,
        "simlen": 200,
    }
    params = fwdpy11.ModelParams(**pdict)

    return pop, params

Tracking a mutation for a specified number of generations#

ALPHA = -10.0

rng = fwdpy11.GSLrng(12345)
pop, params = setup()

mutation_data = fwdpy11.conditional_models.NewMutationParameters(
    frequency=fwdpy11.conditional_models.AlleleCount(1),
    data=fwdpy11.NewMutationData(effect_size=ALPHA / 2 / pop.N, dominance=1),
    position=fwdpy11.conditional_models.PositionRange(left=0.49, right=0.51),
)

output = fwdpy11.conditional_models.track_added_mutation(
    rng, 
    pop,
    params,
    mutation_data,
    when=3,
    until=7,
)

When tracking deleterious variants, it is unlikely that they will be around at the end of the simulation:

if output.mutation_index is not None:
    print(output.pop.mutations[output.mutation_index])
elif output.fixation_index is not None:
    print(output.pop.fixations[output.fixation_index])
else:
    print("Our mutation is no longer in the population!") 

Our mutation is no longer in the population!

Recording all generations of the mutation’s sojourn#

So, we’ve lost all the information about this variant. That’s not so useful. Let’s record all generations of its existence as ancient samples:

rng = fwdpy11.GSLrng(12345)

output = fwdpy11.conditional_models.track_added_mutation(
    rng, 
    pop,
    params,
    mutation_data,
    when=3,
    until=7,
    sampling_policy=fwdpy11.conditional_models.AncientSamplePolicy.DURATION,
)

Now, our mutation is present in nodes in our tree sequence. Let’s try to print it again:

if output.mutation_index is not None:
    print(output.pop.mutations[output.mutation_index])
elif output.fixation_index is not None:
    print(output.pop.fixations[output.fixation_index])
else:
    print("Our mutation is no longer in the population!") 

Mutation[position:0.503040, effect size:-0.010000, dominance:1.000000, origin time:3, label:0]

Let’s track this variant’s frequency at each time point:

for time, nodes, _ in output.pop.sample_timepoints(include_alive=False):
    print(time, len(nodes))
    tree_itr = fwdpy11.TreeIterator(output.pop.tables, nodes)
    for tree in tree_itr:
        for mutation in tree.mutations():
            print(time, tree.leaf_counts(mutation.node), mutation.key)

4.0 1000
4.0 2 0
5.0 1000
5.0 1 0
6.0 1000
6.0 1 0

7.0 1000
7.0 1 0