Calculating number of extant lineages over time

[JaneeshBansal]

Hi,

I have been running simulations in msprime and want to analyse my tree.

For example: If I have a simple simulation ...

ts = msprime.sim_ancestry(1, recombination_rate=0.9, sequence_length=10, random_seed=2)
SVG(ts.draw_svg())

Which produces this tree ...

Is there already a function in tskit that will output the number of extant lineages over time? Or if I could know when the number of branches changes and by how much?

What I thought ...

• I could use nodes_time to get all the times when a coalescent event occurs. Subset to remove the 0's. Then at each node time I could subtract 1 from the total number of lineages I had at the start (in this case 16) and continue until I reach 8. However, this does not work if I have the same node time in two segments (ID 6) as only the unique nodes_time are shown. Is there a way to get all the nodes time?

I may be overcomplicating this, but I would like to calculate the coalescent rate over time for a particular simulation and I think this would help me do so.

Thank you!
Janeesh

[jeromekelleher]

Getting the number of extant lineages overall is a bit tricky, I'm not sure we have an easy way to do that. Related methods which might be helpful:

• Tree.num_lineages gets the number of lingeages in a tree at a given time
• decapitate cuts off the trees sequence at a given time.

This might be an example of an question that can be answered more easily using a graph-based approach like discussed here: #2869

@GertjanBisschop, do you know of an easy way to do this?

[GertjanBisschop]

Hmm, not sure if I myself would think of this particular starting situation as representing 16 lineages. I would probably argue there are two until you hit the first recombination event.
I guess if you really want to count the number of branches in each local tree as observed through time you could use the edge_diffs() (see docs here). You could start out by initializing an array that is the size of np.unique(ts.nodes_time).size which will keep track of the number of branches in tree i. The edge_diffs will give you the branches that are removed to go from one local tree to the next. For each of those you decrement the corresponding interval. Similarly, for each incoming edge you increment the corresponding interval. Before moving to the next tree you then add this intermediate array to an array of the same size that keeps track of the total.

I can point you to some code that does something very similar to what I described above. But I might be completely misunderstanding what you want to achieve.

[JaneeshBansal]

Thank you for the replies!

Would that method give the time at which the number of branches has changed and by how much? I thought edge_diffs() only provides which branches are removed or added into the adjacent local tree compared to its neighbour.

I would like to get the number of branches over the entire sequence (all the trees) over the generations, so the frequency of branches over time. Or at what time the number of branches changes and by how much it changes (as sometimes there are multiple coalescence events occurring during that particular time). I was trying to keep all of the nodes_time even the duplicates to store how many coalescence events were occurring but nodes_time() only keeps the unique times besides the sampled leaves.

I have drawn what I want to extract out of the tree sequence. So I would like t1, t2... with the number of branches that existed during that period. The output would be nice as a list of tuples.

I saw in msprime that the IICR could be calculated using the demographic model specified, I was wondering if the same idea could be applied to a simulation?

I hope this makes sense

[nspope]

I saw in msprime that the IICR could be calculated using the demographic model specified, I was wondering if the same idea could be applied to a simulation?

Hmm, are you wanting to calculate an empirical pairwise coalescence rate from a simulated tree sequence, within specified time intervals? If so, one way to do so is:

import msprime
import numpy as np

ts = msprime.sim_ancestry(100, recombination_rate=1e-3, sequence_length=1e6, random_seed=2)
time_breaks = np.linspace(0.0, 4.0, 11)
coaltime = ts.coalescence_time_distribution(weight_func='pair_coalescence_events', span_normalise=True)
empirical_rate = coaltime.coalescence_rate_in_intervals(time_breaks).squeeze()

print(np.round(empirical_rate, 2)) #theoretical rate is 0.5
# [0.5  0.48 0.56 0.56 0.51 0.41 0.54 0.52 0.54 0.49]

The machinery here is still WIP (undocumented), but the basic idea is to calculate the empirical CDF of pair coalescence events, and use Kaplan-Meier-like estimators for the rate given the ECDF. In the case where the simulation has a single diploid sample, this is effectively the same as estimating the CDF from the ages of internal nodes (e.g. the sole ancestor in each tree), weighting "observations" by the span of each tree. For the case where trees have more than two tips, one can weight nodes (in a given tree) by the number of pairs that coalesce in each node.

This isn't quite the same as calculating the number of extant lineages (although it might be close in the case of two-tip trees) -- but I think a similar approach could be used to do so.

[JaneeshBansal]

Thank you for the code. I have got it running now! I just had to ensure that my window sizes were large enough otherwise enough coalescence events were not occurring.

Janeesh

[nspope]

Yes, the rates are independent estimates per window, akin to interval-censored hazard rate estimators in time-to-event analysis. So, there'll be substantial variance if the windows are too small. One can improve things substantially with small windows by smoothing, but if you're using these on simulations, I think it probably is easiest just to simulation longer sequences / more samples to decrease the variance.

[JaneeshBansal]

Hi,

I was wondering where the code for the functions coalescence_time_distribution() and coalescence_rate_in_intervals() is? As it is undocumented in the msprime manual, where would it be in the GitHub source code?

I think the question I am curious about is how is the coalescence rate is calculated within the specific time breaks selected. For example for generations [0:3000] is an average calculated over this interval for the rate?

Thank you!
Janeesh

[GertjanBisschop]

Hi @JaneeshBansal,

It was lovely meeting you in person at Probgen. Thanks for reminding us we forgot to answer your question. This is the bit of code that does all that estimates the empirical coalescence rate. Perhaps also look at the original pull request If you need further info after reading this, then we'll ping @nspope again.