Create a new Tree() by an SPR from another tree

[hyanwong]

It might be useful to have a method that produces a new tree by a single SPR move from an existing tree. Perhaps the way to do this is to take an entire tree sequence and return a new tree sequence which has been extended in length by (say) 1 unit, with the new tree tacked onto the end (and the edge and node tables extended). That way it would be possible to build up a tree sequence by repeatedly calling this method. It would be sensible to have an equivalent table method too.

It does, however, require some model of the probabilities of generating different SPRs. So perhaps this is too "model based" for a tskit method, and it belongs in something like msprime (and on reflection, perhaps this discussion belongs in the msprime repo instead). I think there is a standard SMC way to make a new tree from an old one, assuming recombination is equally likely on any branch (proportional to branch length).

For context, this arose during a really useful chat to @a-ignatieva on Friday, where we discussed an ARG distance metric which could take account of the persistence of edges between trees, so that a "proper" ARG would have lower (better) score than a JBOT of the same marginal topologies. For testing purposes, it would be useful to create trees that were "marginally" different from existing trees.

[jeromekelleher]

It's a nice idea, but it's not obvious to me what the proposed API would look like. Any suggestions?

[hyanwong]

I guess either tables.append_spr(sequence_length=1) or msprime.append_spr(tables, sequence_length=1), which would increase the sequence length by the given amount and plonk a new tree there?

(we might need other parameters if it turns out, as I suspect, that we need to parameterise the model used to generate the SPR. For example, if the SPR creates a new MRCA for the tree, what time should the new MRCA node have?)

[jeromekelleher]

It depends - are you talking about having a stochastic process in which you add on an extra tree (i.e. implementing the Wiuf and Hein algorithm or the SMC), or are you talking about the mechanistic details of just updating the tables with (essentially) an element of the edges_diffs list?

If it's implementing the stochastic process, this isn't trivial and would need substantial validation.

[hyanwong]

Hmm, good point. On reflection I think there are 3 separate functions here:

1. Take an edge_diff (an set of edges in and edges out) together with a list of any new nodes and append a new tree onto the end of the tree sequence. This needn't correspond to a single SPR.
2. Use 1. to append a new tree that is 1 SPR different from the last one in the tree sequence, by specifying a single node to prune and an insertion point (with associated time) where it can be validly regrafted.
3. Use 2. plus some arguments about branch lengths, models etc. and a pRNG to implement a stochastic SPR addition algorithm.

I think 1 and (probably) 2 could belong in tskit. 3 might simply be a recipe in the tutorials, if the algorithm is easy. Or if could be in a separate library.

Does that make any sense?

[jeromekelleher]

Makes sense, yeah. Would be best to sketch out the full process before trying to add anything to tskit or msprime though I think, there will be tricky details rearing up in 3.

[benjeffery]

This sounds great, but also sounds like a separate repo from both msprime and tskit?

[hyanwong]

Hm, possibly. What would we call it? tsmanip perhaps? I guess the same argument could be said of e.g. the stats functions. I think of it as a tskit thing because I think it's nice to be able to say tables.append_spr() rather than my_lib.append_spr(tables) or new_ts = my_lib.append_spr(ts). But whatever others prefer, really.

[benjeffery]

For experimental things like this, we can always bring them into tskit when they are polished. Makes sense to me to do the exploration and initial dev in a new repo - doesn't have to have all the structure of an official tskit one.