Convert forward simulation notebooks.

[jeromekelleher]

There's two notebooks (wfforward.ipynb and wfcython.ipynb) in the old-content/notebooks directory that need to be converted to jupyterbook format. Hopefully this can be reasonably automatic using jupytext.

The only concern is the environment requirements - what is need to run these notebooks? It would be great if we keep this a no-conda affair.

@molpopgen, would you mind taking a look please?

[jeromekelleher]

Also, we want to make sure these pages don't take to long to run.

[molpopgen]

I doubt they need conda--apt should be fine. The one with cython took quite a while to run, in order to demonstrate correctness. I have previously suggested that we delete it.

[jeromekelleher]

Your call on deleting stuff. Demonstrating correctness is beyond the scope of these tutorials now I think.

[molpopgen]

Looking at it now. It was written back when the access to tskit was via msprime. It'll take some time to convert. Also, there's an annoyance--if you convert with jupytext, you can get YAML frontmatter that somehow is incompatible with what is in other md files, causing the whole thing to fail to build.

[jeromekelleher]

I've not bothered modernising the content in the first pass @molpopgen, if that's any help. I've basically just been getting things working, and going to make notes on the outstanding issues so we can follow up in smaller chunks afterwards.

[molpopgen]

yeah, I'm going for "working" right now, which it doesn't!

[molpopgen]

I'm going to shorten it considerably. Looking at it now, I really don't like the length.

[molpopgen]

This will take a bit more time than I thought. Teaching, etc.,..

[hyanwong]

Hi @molpopgen - I'm just looking at our missing tutorials, and see there's a space for forward simulations. Did you get anywhere with this?

[molpopgen]

No, totally off of my radar. I still vote to delete the Cython example and then minimally port over the other to markdown.

[hyanwong]

No, totally off of my radar. I still vote to delete the Cython example and then minimally port over the other to markdown.

That sounds ideal to me. You mentioned something about shortening though? Is that still relevant?

[molpopgen]

That sounds ideal to me. You mentioned something about shortening though? Is that still relevant?

Yes, IMO a tutorial should strive to be only a "screen or 2 long". To that end, something like cell 3 modified to do exactly 1 crossover / meiosis in the "pseudo-haploid" style of the C example will be sufficient. Updating to include periodic simplification would be a good thing, too, I think.

[hyanwong]

IMO a tutorial should strive to be only a "screen or 2 long"

@molpopgen - I am looking at this with a view to showing a simple forwards simulator for a workshop. I rethought the python code from scratch, and came up with this. I would be interested in how clear or otherwise you think it is (I am guessing it isn't very efficient though). The nice thing about picking parents, rather than genomes, at random is not just that it is more obvious to a biologist, but that we can relatively easily modify the add_inheritance_paths function below to allow recombination between the two genomes within a parent (EDITED to make nicer code):

import tskit
import random

# First building block: make an individual containing 2 genomes. Note that
# we flag up the nodes as "sample nodes", which helps when plotting.
def make_individual(tables, time):
    """
    Creates a new individual containing 2 new genomes (nodes), and returns the ID
    of the individual and a tuple of the IDs of the individual's two genomes 
    """
    indivID = tables.individuals.add_row()
    mat_nodeID = tables.nodes.add_row(tskit.NODE_IS_SAMPLE, time, individual=indivID)
    pat_nodeID = tables.nodes.add_row(tskit.NODE_IS_SAMPLE, time, individual=indivID)
    return indivID, (mat_nodeID, pat_nodeID)

# Second building block: add "inheritance paths" between the
# genome the child inherits from one parent (say the mother) and
# the two genomes contained in that parent. This is where "meiosis"
# comes in. However the moment, we ignore recombination and simply pick
# one of the two genomes at random

def add_inheritance_paths(
    tables, parent_genomes, child_genome, recombination_rate=None
):
    """
    Given two parent genome (node) IDs from the same individual, add inheritance
    paths from those genomes to a child genome (node) ID. 

    """
    left = 0
    right = tables.sequence_length
    inherit_from = random.randint(0, 1)  # randomly chose 0 or 1
    tables.edges.add_row(left, right, parent_genomes[inherit_from], child_genome)

# Last building block: take a collection of parents 
# and generate children from them by picking two parents at random
# Here we use `random.choices` to sample with replacement, but we could
# use `random.sample` to exclude the possibility of selfing
def make_new_generation(tables, time, num_children, parents):
    """
    Returns a Python dictionary of children (i.e. length `num_children`) mapping
    {individual_ID: (maternal_genome_ID, paternal_genome_ID), ...}
    The `parents` parameter is a dictionary of the same form, listing possible parents
    for the children. If the dictionary is empty, no inheritance paths are set.
    """
    children = {}
    parent_IDs = list(parents.keys())
    for _ in range(num_children):
        child_ID, child_genomes = make_individual(tables, time)
        children[child_ID] = child_genomes
        if parents:
            # pick IDs for two parents at random with replacement (allows selfing)
            random_parent_IDs = random.choices(parent_IDs, k=2)
            for child_genome, parent_ID in zip(child_genomes, random_parent_IDs):
                add_inheritance_paths(tables, parents[parent_ID], child_genome)
    return children

def setup_simulation(sequence_length, random_seed):
    "Returns a TableCollection where we can store our growing genealogy"
    random.seed(random_seed)
    tables = tskit.TableCollection(sequence_length)
    tables.time_units = "generations"
    return tables

def simple_diploid_sim(popsize, generations, sequence_length, random_seed=123):
    tables = setup_simulation(sequence_length, random_seed)

    # tskit time is in units ago: we want the final generation to be at time 0
    generations_ago = generations - 1  # start with the oldest generation ago
    
    # Initialise an empty population of parents
    parents = make_new_generation(tables, generations_ago, popsize, parents={})

    # Run the simulation
    while generations_ago > 0:
        generations_ago = generations_ago - 1
        parents = make_new_generation(tables, generations_ago, popsize, parents)

    # Sort the tables, e.g. to put edges in the required order
    tables.sort()
    return tables.tree_sequence()

ts = simple_diploid_sim(5, generations=15, sequence_length=1000)

Here's how you might swap in a function that allows recombination:

def add_inheritance_paths(tables, parent_nodes, child_node, recombination_rate=0):
    """
    Add inheritance paths from parent node IDs (i.e. parent genomes)
    to a child node ID (i.e. a child genome)
    """
    breakpoints = ...  # e.g. pick breakpoints with probability rho per integer position
    inherit_from = random.randint(0, 1)
    for lft, rgt in zip([0] + breakpoints, breakpoints + [tables.sequence_length]):
        tables.edges.add_row(
            left=lft, right=rgt, parent=parent_nodes[inherit_from], child=child_node)
        inherit_from = 1 - inherit_from  # switch: inherit from the other parent genome

[molpopgen]

@hyanwong I find it a bit hard to read for a few reasons:

• variable names differ in a few places from function argument names. parents vs parent_genomes for example. For a "tutorial", I think it is best for names to be consistent throughout.
• I'd suggest type hints throughout as well.
• make_new_generation is obfuscated. Passing in the empty dictionary the first time is non-obvious. I think the input type would be Optional[Dict[...]] and pass in None the first time.