Checking whether a node has mutations above it in a given tree

[a-ignatieva]

The premise is that I go to a tree T1 at a specific site, select a (small) group of nodes from this tree based on some metric, and then move to tree T2 just to the right of T1. For each node in my group, I would like to check whether the node has any mutations (within the span of T2) on the edge just above it. Is there a quick and efficient way of doing this?

[hyanwong]

One way to do this is to iterate over all the mutations in the new tree (very efficient), and check if the mutation.node is in your set. Or you could grab all the mutation node IDs as a numpy array and do it at ones using np.isin. That should be fast as long as there are not bazillions of mutations per tree, I think?

[hyanwong]

The mutations should be ordered by site ID, which itself is ordered by position, so all the mutations in a tree will be contiguous in the mutations table. That means you should be able to take a slice of ts.tables.mutations.nodes to grab all nodes with a mutation above them in a given interval. However, finding the left and right indexes of the slice is a bit involved, and it's going to be a lot simpler to do np.array([m.node for m in tree.mutations()]). So I'm not sure it's worth the bother of working out the slice positions. Perhaps others will have an opinion?

[jeromekelleher]

I would do it the simple way first as Yan suggests above. Let S be the set of nodes you're interested in, then in T2 do

for site in T2.sites:
    for mutation in site.mutations:
          if mutation.node in S:
              print("Mutation in T2 at node {mutation.node}")

If you're doing this a lot (and you have lots of trees along the genome) then how you iterate over the trees can be important. Happy to provide example code here if you'd like to be more specific @a-ignatieva

[a-ignatieva]

Thanks @hyanwong and @jeromekelleher!

Let me try this out and compare it to slicing the table as Yan suggests - my trees do tend to have lots of mutations, so maybe the latter could be faster? I have something like the following for finding the left and right sites that correspond to a given genomic interval, which I could use to find the left/right indices of the slice I need - though not sure if this is the best way of doing this?

def pos_to_site(ts, interval):
    if interval[0] < ts.tables.sites.position[0]:
        a1 = -1
    elif interval[0] > ts.tables.sites.position[-1]:
        a1 = ts.num_sites
    else:
        a1 = np.min(np.where(ts.tables.sites.position >= interval[0]))
        
    if interval[1] < ts.tables.sites.position[0]:
        a2 = -1
    elif interval[1] > ts.tables.sites.position[-1]:
        a2 = ts.num_sites
    else:
        a2 = np.max(np.where(ts.tables.sites.position <= interval[1]))
        
    return [a1, a2]

I also do iterate over a lot of trees - I basically keep moving one tree to the right until some condition is reached. I use T.next() for this - is there a more efficient way?

[hyanwong]

That gets you the IDs of the first and last site in the tree, but you want the first and last mutation in the tree, which would be a little bit more involved, I'm afraid. Also your np.min(np.where(ts.tables.sites.position >= interval[0])) is looking through the entire sites table, so there's probably a way to do this so that it doesn't search through the entire thing in each loop.

But I suspect that this may be premature optimization. The method @jeromekelleher pointed to above should be very fast, especially if you can wrap it in a numba decorator (I've not tried this, see tskit-dev/tutorials#63)

[jeromekelleher]

I agree with Yan - I'd try it the simple way first and do some profiling if it's too slow. If you have some code we could run on some simulated data I'd be happy to take a look.

The method @jeromekelleher pointed to above should be very fast, especially if you can wrap it in a numba decorator (I've not tried this, see tskit-dev/tutorials#63)

Numba won't help here because we're using Python APIs, you'd need to be working directly with numpy arrays for it to make a big difference.

I also do iterate over a lot of trees - I basically keep moving one tree to the right until some condition is reached. I use T.next() for this - is there a more efficient way?

No, that's as good as you're going to do @a-ignatieva!