How to find the total genome contributed by each founder?

[petrelharp]

Over on the SLiM list, Savanah Bird asks:

I'm using SLiM to simulate a wildlife population re-introduction. I'm burning in a source population, then selecting a handful of individuals to be the founders of a new population, then letting the simulation run forward a short amount of time (25 ticks). In my simulations, I'd like to calculate the founder genome equivalents (Lacy 1989) of the re-introduced population, which takes into account 1) the proportion of each founder's genome that has been retained by the extant population (i.e. not lost due to drift) and 2) the proportion of total genomic material in the extant population contributed by each founder.

Here's some code that does that, by iterating over the roots in each tree. For this script to work we need to Remember each of the founder individuals (else it's not as easy to find which chromosomes are in the same individual).

import tskit
import pyslim
import numpy as np

ts = tskit.load("founders.trees")

def founder_stats(ts, root_time):
    """
    Returns arrays with (a) the total amount of genome on which
    each founder has left genetic material, and (b) the total amount
    of genetic material in the final population. This is calculated
    as a proportion of the genome, averaged across chromosomes,
    so a diploid founder who has descendants carrying one chromosome
    but not the other will have "has_offspring" of 0.5.
    """

    founder_indivs = pyslim.individuals_alive_at(ts, root_time)
    modern_indivs = pyslim.individuals_alive_at(ts, 0.0)
    modern_nodes = np.array([n for i in modern_indivs for n in ts.individual(i).nodes])

    node_indexes = np.full((ts.num_nodes,), -1, dtype='int')
    for k, i in enumerate(founder_indivs):
        ind = ts.individual(i)
        for n in ind.nodes:
            node_indexes[n] = k

    has_offspring = np.zeros(len(founder_indivs))
    total_offspring = np.zeros(len(founder_indivs))
    for t in ts.trees(tracked_samples=modern_nodes):
        for r in t.roots:
            k = node_indexes[r]
            assert k >= 0, f"{r} not in list"
            f = t.num_tracked_samples(r)
            print(k, f, t.span)
            has_offspring[k] += t.span * (f > 0)
            total_offspring[k] += t.span * f

    has_offspring /= 2 * ts.sequence_length
    total_offspring /= 2 * ts.sequence_length
    return has_offspring, total_offspring

root_time = np.max(ts.tables.nodes.time)
has_offspring, total_offspring = founder_stats(ts, root_time)

Here's the example SLiM script to do this:

initialize() {
	initializeTreeSeq();
	initializeMutationRate(0.0);
	initializeMutationType("m1", 0.5, "f", 0.0);
	initializeGenomicElementType("g1", m1, 1.0);
	initializeGenomicElement(g1, 0, 9999999);
	initializeRecombinationRate(1e-8);
}

1 early() {
	sim.addSubpop("p1", 4);
	sim.treeSeqRememberIndividuals(p1.individuals);
}
5 late() {
	sim.treeSeqOutput("founders.trees");
}

[savanahbird]

Hi Peter,

I just tested this out, and it works! The function is just missing a return statement (it'll give a "cannot unpack non-iterable NoneType object" as is).

Thanks again!

[petrelharp]

thanks! fixed the return statement.

[savanahbird]

Building off founder_stats function to calculate Founder Genome Equivalents:

def calc_fge(ts, root_time):
    """
    For a given tree sequence and root time, returns the founder genome equivalents
    as defined by Lacy 1989. Calls founder_stats function, which returns arrays with 
    (a) the total amount of genome on which each founder has left genetic material, 
    and (b) the total amount of genetic material in the final population.
    """
    has_offspring, total_offspring = founder_stats(ts, root_time)
    terms=[]
    for i in range(len(has_offspring)):
        r=has_offspring[i]
        if r>0:
            p=total_offspring[i]/sum(total_offspring)
            terms.append(p**2/r)
    return float(1/sum(terms))
        

fge=calc_fge(ts, root_time)