Eliminate only unary nodes beyond the TMRCA of each tree

[RomainPechey]

Hello, I am trying to count the number of ancestors separating each pair of individuals i,j. To do this, I ran simulations in SLiM and extracted tree sequences containing unary nodes.
I would now like to represent these trees, but I am unable to remove the unary nodes located beyond the TMRCA of each tree (the “tail” of each tree).
I looked at the arguments of the simplify() function.
When I use keep_unary=True, all unary nodes are retained; when I switch to keep_unary=False, the trees are cut correctly, but I lose the unary nodes on my branches.
I also tried keep_input_roots=False, but that didn't change anything. Is there a simple function that would solve my problem?
I am attaching an example of a tree that I would like to simplify.
Thank you in advance.

Romain Pechey

[hyanwong]

Yes, the keep_unary option is rather crude at the moment. One way you could do this is to pass all the nodes you want to keep (including all the unary ones), but use the update_sample_flags=False option:

import msprime
import numpy as np

ts = msprime.sim_ancestry(50000, sequence_length=1000, recombination_rate=0.0002, random_seed=1234).simplify(np.arange(4), keep_unary=True)
ts.draw_svg(time_scale="log_time")

drop = np.ones(ts.num_nodes, dtype=bool)
node_is_sample = np.zeros(ts.num_nodes, dtype=bool)
node_is_sample[ts.samples()] == True
for tree in ts.trees():
    # find the first node below the root that has 2 children
    u = tree.root
    num_children = tree.num_children_array
    left_child = tree.left_child_array
    while num_children[u] == 1 and not node_is_sample[u]:
        u = left_child[u]
    drop[tree.preorder(u)] = False
new_ts = ts.simplify(np.where(drop == False)[0], update_sample_flags=False)
new_ts.draw_svg()

[hyanwong]

However, note that this keeps e.g. node 80 in the 2nd tree, because it is a root in the first tree. You might not want this?

[RomainPechey]

Thank you very much for your reply. Unfortunately, it does not fully address my issue. In fact, I am trying to do this in order to determine the “actual” number of marginal trees contained in my trees. However, after running a few tests, I found that unsimplified tree sequences contain more marginal trees than simplified trees (properly cut at the root).Upon further investigation, it seems to me that some trees in the unsimplified tree sequence share the same topology, but only the “tail” part differs. In other words, I feel that this tail “artificially” increases the number of marginal trees in the tree sequence.Am I mistaken? (I am attaching another tree sequence as an example: trees 2 and 3 have exactly the same topology, only the tail differs, as do the penultimate and last trees.).

P.S : I am one of the doctoral students working with Pierre-Alexandre Gagnaire on the use of IBD for demographic inference. I would like to refer you to the code that we have posted online with Iago Bonnici (@iago-lito) and Océane Eychenne (@OceaneEYCHENNE) [(https://github.com//discussions/3142),] which currently allows us to obtain all the stitched IBD pairs, while we wait for this step to be bypassed.

[hyanwong]

Thanks, and good to know the background. It could easily be that there are recombination events above the local MRCAs, and these are creating more trees. If you only want to count the number of trees in the tree sequence, why don't you just do that after full simplification? You can always go back to the unsimplified version if you need to count the number of intermediate nodes (or just use time, right?)

If you are interested in IBD, you might want to simplify but keep the unary portions of coalescent nodes (e.g. see #2127 which has some code to hack this). Alternatively you could fully simplify, then extend_haplotypes, but that will give a slightly different answer (and probably fewer trees).