Function to return the common ancestor of a set of nodes

[stsmall]

Hi all,
Is there a function in tskit to return the lowest common ancestor of a set of nodes? For example, the lowest common ancestor of a population. I tried writing something using networkx function all_pairs_lowest_common_ancestor, where I keep iterating through the returned mrca until only 1 remains … but it is pretty inefficient. Any ideas?
thanks!

[jeromekelleher]

There isn't a function to do it currently in tskit @stsmall, and I guess we probably should add it.

In the mean time, this will work:

import tskit
import functools

tree = tskit.Tree.generate_balanced(9, arity=2)
print(tree.draw_text())

nodes = [0, 2, 3]
mrca = functools.reduce(tree.mrca, nodes)

print(mrca)

gives:

       16        
   ┏━━━━┻━━━┓    
   ┃       15    
   ┃     ┏━━┻━┓  
  11     ┃   14  
 ┏━┻━┓   ┃  ┏━┻┓ 
 9  10  12  ┃ 13 
┏┻┓ ┏┻┓ ┏┻┓ ┃ ┏┻┓
0 1 2 3 4 5 6 7 8

11

This works because we can successively apply Tree.mrca(nodes[0], Tree.mrca(nodes[1], nodes[2])).

We should probably add this to the documentation, since it's such a neat and simple recipe.

[jeromekelleher]

Asymptotically, I think this is probably just as efficient as tracking the paths through the tree at the same time, but in practise it won't be as fast because of the Python overhead. We could easily extend the Tree.mrca function to take more positional arguments, and then pass the entire list to the C function. I guess we could do a pairwise reduction in C without needing to do anything fancy then, algorithm-wise.

[jeromekelleher]

If someone was doing a very large number of MRCA calculations, then it would probably be worthwhile doing it in numba using the magical Schieber-Vishkin algorithm (implementation here)

[petrelharp]

@gtsambos had a good suggestion also - to use simplify(). If you need the MRCA of the same set of nodes over the entire tree sequence, this might be more efficient than doing it tree-by-tree? So, this would work something like this:

ts = msprime.simulate(10, recombination_rate=0.2, random_seed=7)
sts, nm = ts.simplify([0,4,6], map_nodes=True)
inverse_map = np.repeat(-1, sts.num_nodes)
for j, k in enumerate(nm):
    if k != tskit.NULL:
        inverse_map[k] = j

sts_trees = sts.trees()
st = next(sts_trees)
for t in ts.trees():
   if st.interval[1] <= t.interval[0]:
     st = next(sts_trees)
   print(t.draw(format='unicode'))
   print(f"The MRCA is {inverse_map[st.root]}.")

which produces

       18          
    ┏━━━┻━━━┓      
   17       ┃      
 ┏━━┻━┓     ┃      
 ┃    ┃    16      
 ┃    ┃   ┏━┻━━┓   
15    ┃   ┃    ┃   
┏┻┓   ┃   ┃    ┃   
┃ ┃  14   ┃    ┃   
┃ ┃  ┏┻━┓ ┃    ┃   
┃ ┃  ┃  ┃ ┃   13   
┃ ┃  ┃  ┃ ┃  ┏━┻━┓ 
┃ ┃  ┃  ┃ ┃  ┃  12 
┃ ┃  ┃  ┃ ┃  ┃  ┏┻┓
┃ ┃ 11  ┃ ┃  ┃  ┃ ┃
┃ ┃ ┏┻┓ ┃ ┃  ┃  ┃ ┃
┃ ┃ ┃ ┃ ┃ ┃ 10  ┃ ┃
┃ ┃ ┃ ┃ ┃ ┃ ┏┻┓ ┃ ┃
0 8 3 6 4 1 2 7 5 9

The MRCA is 17.
     19            
 ┏━━━━┻━━━━┓       
 ┃        18       
 ┃    ┏━━━━┻━━━━┓  
 ┃   16         ┃  
 ┃  ┏━┻━━┓      ┃  
15  ┃    ┃      ┃  
┏┻┓ ┃    ┃      ┃  
┃ ┃ ┃    ┃     14  
┃ ┃ ┃    ┃     ┏┻━┓
┃ ┃ ┃   13     ┃  ┃
┃ ┃ ┃  ┏━┻━┓   ┃  ┃
┃ ┃ ┃  ┃  12   ┃  ┃
┃ ┃ ┃  ┃  ┏┻┓  ┃  ┃
┃ ┃ ┃  ┃  ┃ ┃ 11  ┃
┃ ┃ ┃  ┃  ┃ ┃ ┏┻┓ ┃
┃ ┃ ┃ 10  ┃ ┃ ┃ ┃ ┃
┃ ┃ ┃ ┏┻┓ ┃ ┃ ┃ ┃ ┃
0 8 1 2 7 5 9 3 6 4

The MRCA is 19.

[stsmall]

Why does one have to use map_nodes in this case?

[petrelharp]

That's how we know how to associate the nodes in the simplified tree sequence with those in the original ones. If all you want is some attribute of theirs (e.g., their time), you don't need to do that.

[benjeffery]

Thanks for the question @stsmall - we've added an issue #1340 that will bring this functionality into tskit!

If you're motivated to do so you could author the pull request yourself and get credit as a contributor. Guides to do so are at https://tskit.dev/tskit/docs/stable/development.html . You'll only need to touch the python code so ignore sections of that document that refer to the C API.

[stsmall]

Hi All,
Thanks for the great suggestions!
A related question ... I was reading Hejase et al. 2020, and they use a statistic, TMRCAH, which was defined as the smallest (shortest) subtree that contains half the nodes of a population. What would be an efficient way to calculate the TMRCAH statistic? Would this be a case where I would need to trace all the paths from the leaves to the root ... as Daniel Goldstein suggested in the slack discussion? Maybe combining that with Peter/Georgia suggestion of using simplify?

[petrelharp]

This seems like a good topic for a different "discussion" , to keep things easily findable? And, a very good question...

[stsmall]

OK, np. I will open it. Thanks!

[hyanwong]

This can now be done by https://tskit.dev/tskit/docs/stable/python-api.html#tskit.Tree.mrca. This should probably be selected as the "correct" answer now.