Tracking user-specified new mutations

Show code cell source

Hide code cell source

import fwdpy11
import numpy as np
import msprime

import fwdpy11.conditional_models

Show code cell source

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,
        "demography": fwdpy11.ForwardDemesGraph.tubes([pop.N], burnin=10)
    }
    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