get shared tree across an interval

[petrelharp]

A common thing we'd like to do is to visualize shared structure in the trees across a region of the genome. It occurs to me that one way to do this is would be to be able to do like ts.at( ), perhaps

tree = ts.at(interval=(left, right))

which returns the tree constructed by all edges that span the entire interval [left, right).

Probably, this would be like "if the argument to ts.at( ) is a tuple of length two, return the interval".

I do like the ts.at() syntax, but we might alternatively implement this as

def tree(at=None, interval=None):
    if at is not None:
        return self.at(at)

and deprecate ts.at().

[hyanwong]

Isn't this the same as having a left and right parameter for trees()?

There's a PR I made to allow this: #2762

[petrelharp]

Isn't this the same as having a left and right parameter for trees()?

I don't think so? The way I read #2762 is that it gives you a way to iterate over all the trees from left to right. What I'm proposing here is a way to return the single tree that consists of all edges that are shared among all those trees between left and right. So, this resulting tree will have multiple roots unless there's only a single tree that spans left to right.

Here's an implementation:

import msprime, tskit
import numpy as np

ts = msprime.sim_ancestry(2, recombination_rate=1e-8, sequence_length=1e8, random_seed=123)

def tree(ts, left, right):
    tables = ts.dump_tables()
    tables.edges.keep_rows(
            np.logical_and(
                tables.edges.left <= left,
                tables.edges.right >= right
            )
    )
    nts = tables.tree_sequence()
    # this is only for error checking
    if nts.num_edges > 0:
        assert nts.keep_intervals([[left, right]]).trim().num_trees == 1
    return nts.at(left)

which gets

# print(tree(ts, 0, 10).draw(format='unicode'))
 
  11   
 ┏━┻━┓ 
 9   ┃ 
┏┻┓  ┃ 
┃ ┃  4 
┃ ┃ ┏┻┓
0 3 1 2

# print(tree(ts, 0, 1e7).draw(format='unicode'))
  11   
   ┃   
   4   
  ┏┻┓  
0 1 2 3

[jeromekelleher]

Ahh, now I get it. This is cool - it needs a better name though. Ideas

• ts.shared_tree(left, right)
• ts.spanning_tree(left, right) (but confusing for graph theory people, unless it's a good analogy?)
• ts.constant_tree(left, right) (but a bit misleading because abutting edges going in and out will not be included)
• Something else? I don't think overloading ts.at is the right way to go anyway though, it's quite a distinct thing IMO

[hyanwong]

How about ts.shared_edges or ts.shared_edge_tree?

[petrelharp]

I like ts.shared_tree (assuming this returns a Tree)!