all nodes_time() - even duplicates

[JaneeshBansal]

Hi,

I was wondering if it is possible to get all of the nodes_time() from a tree? As currently it removes the duplicates - for example, could I get the node time for ID 4 written out three times?

I would like to produce a histogram of the coalescence time distributions for my simulation.

Thank you!
Janeesh

[petrelharp]

Hm - to get the distribution of coalescence times wouldn't you want to weight be sequence length, not number of trees?

If so, @nspope has written some methods to do this.

[petrelharp]

Well, the short answer to your question is that you could do like:

times = []
for t in ts.trees():
  for n in t.nodes:
    times.append(ts.node(n).time)

... but I'm not sure that's what you actually want to do?

I think one question here is what you want to mean by the "frequency of coalescence event", and I know of two natural interpretations:

1. One event is one distinct coalescence between ancestral segments, and you want the frequency of these. In this case, the same event is in general shared across many trees (and more recent events across more trees). This would be something like the histogram of node times (but this is not right in general).
2. Or, the distribution of coalescence time, averaged across locations on the genome - this is what is predicted by coalescent theory, for instance - and for this, you'd want to weight by sequence length. This is what the currently undocumented CoalescenceTimeDistribution does - but it's not ready for general use.

So - in neither of these interpretations do you want to weight by the number of trees they appear in. But, if that's actually what you want, then great!

[nspope]

As Peter says, depends on what you want. If you want to estimate the "coalescence rate" or "density of coalescence times" of the marginal process -- marginal in the sense that a random draw from the process could be obtained by picking a genomic coordinate at random and extracting the genealogy at that coordinate -- then you'll need to weight each tree by the probability of sampling it (which is the relative span).

[JaneeshBansal]

Sorry I am still developing my understanding of coalescent theory so these may be silly questions!

I think my initial confusion was about those individual segments parent being ID4. I thought that they were all independent coalescence events, but after reading the msprime paper I think those segments would have a common ancestor which is the union of those three segments? Then that common ancestor would undertake a coalescent event. But in my mind, I thought I had 3 individual coalescent events so was confused as to why they were not being recorded.

I think I would like to know the coalescence times weighted by the length of the sequence. Could you let me know how I can do this please?

[nspope]

No worries, these aren't silly questions at all. If you're interested in pairwise coalescence events, and you're simulating a single diploid, then it's very straightforward: your observations and weights are,

import msprime
ts = msprime.sim_ancestry(samples=1, population_size=10000, sequence_length=10e6, recombination_rate=1e-8)
obs = []
weight = []
for tree in ts.trees():
   obs.append(tree.tmrca(0, 1))
   weight.append(tree.span / ts.sequence_length)

You can check that the (appropriately weighted) mean is close to the haploid population size used in the simulation (20,000 in this case), as expected under the coalescent:

import numpy as np
print(np.sum(np.array(obs) * np.array(weight))) # should give 2*population_size asymptotically
# 20590.54

if you're simulating more than a single diploid, it's a little more complicated (especially to do efficiently) -- but does the above get you what you want?

[nspope]

I think those segments would have a common ancestor which is the union of those three segments?

as an aside, I don't quite understand what you mean here ... what's going on with the diploid simulation is that segments have different TMRCAs (e.g. inherit from different ancestral haplotypes), because of recombination. The length of the segments and the TMRCAs are not independent, which is why one has to weight by segment length to get the correct marginal expectation.