Extracting branch lengths from mutations

[Schraiber]

Hi all,

The short story is that I'm trying to learn something about branch lengths under selection, so I've simulated some data in slim, then kept the tree seq, and I'm importing it into tskit in python. I'm by no means a tskit expert (and honestly I'm a pretty bad python coder in general, since R is more my speed), as I've hardly used it at all, but what I'm doing is essentially this:

edges = trees.tables.edges
nodes = trees.tables.nodes
muts = trees.tables.mutations
sites = trees.tables.sites
branch_length = []
for m in muts:
    n = m.node
    s = m.site
    p = sites[s].position
    edge = edges[(edges.child==n)*(edges.left<=p)*(edges.right>p)]
    branch_length.append(nodes[edge.parent].time-nodes[edge.child].time)

which may not be optimal, but runs fast enough for my purposes. My logic here is that the documentation says that the the node member of mutation is the first node "below" the mutation, which I take to mean the node closer in time to the present. Then, I figure out what edges have that node as a child (which again, I take the child of the edge to be the node closer to the present) and also cover the position of the mutation. Then I'm just trying to calculate the branch length by taking the difference in the parent and child times.

The problem is that it doesn't work! I often end up with the edge object being empty, and it seems to be one of two problems:

1. NO edge has n as a child
   or
2. Some edges have n as a child, but none of them span the position of the mutation.

Clearly this means I'm misunderstanding something about how this data is stored, and if you could give me a hand I'd really appreciate it!

I can't seem to upload the .trees file since it's not a supported file type, but happy to share it and the slim code that generated it if necessary.

[nspope]

Here's one approach:

import msprime
import numpy as np
import tskit
ts = msprime.sim_ancestry(samples=100, sequence_length=1e6, population_size=1e4, recombination_rate=1e-8)
ts = msprime.sim_mutations(ts, rate=1.25e-8)
branch_lengths = np.full(ts.num_mutations, np.nan)
for m in ts.mutations():
    if m.edge != tskit.NULL: # skip unmapped mutations; these'll be nan in "branch_lengths"
        p, c = ts.edges_parent[m.edge], ts.edges_child[m.edge]
        branch_lengths[m.id] = ts.nodes_time[p] - ts.nodes_time[c]

[Schraiber]

This is much cleaner and nicer code, thanks! I see that the edge member of the mutation object is definitely the way to go compared to the table-based way I was doing things.

It looks like the things I'm missing don't have edge as nan, but instead have edge as -1---I'm guessing that this is due to the genealogy not coalescing over the course of the slim simulation, which explains why the node was not the child of any edge. I figured I was running the slim sims for long enough, but I guess not! I'll just increase that and see if it solves my problem, which I suspect it will if I'm understanding things correctly.

[nspope]

Hm, off the top of my head I don't know why you'd have mutations where edge is tskit.NULL from a SLiM simulation -- but maybe @petrelharp knows?

[Schraiber]

Yeah, actually looking a little deeper, I'm actually not sure I understand what's going on, for example this mutation

Mutation(id=0, site=0, node=17928, derived_state='580483', parent=-1, metadata={'mutation_list': [{'mutation_type': 1, 'selection_coeff': 0.0, 'subpopulation': 1, 'slim_time': 28966, 'nucleotide': -1}]}, time=11034.0, edge=-1)

is actually pretty late in my slim simulation, which has a pop size of 1k and I ran the sim for 40k generations (and this mutation appeared in generation 28966), so it seems hard to believe it hadn't coalesced, but it has a null parent and a null edge. Again, it's very possible I'm not understanding something about the tskit representation?

Edit: just realized that parent is relevant for stacked/recurrent mutations, not coalescence, so that doesn't matter, but having edge = -1 still seems wrong?

[nspope]

You could look at the tree at site 0 and see where mutation 0 / node 17928 lie (is the mutation above root or isolated subtree?)

m = 0 # mutation id
s, n = ts.mutation(m).site, ts.mutation(m).node
# don't draw mutations/nodes except for the node/mutation in question (red circle; blue cross)
style = ".node > .sym {fill : none} .mut > .sym {stroke : none}"
style += ".n" + str(n) + " > .sym {fill : red}"
style += ".m" + str(m) + " > .sym {stroke : blue}"
ts.at(ts.sites_position[s]).draw_svg(
  path=f"site{s}_node{n}.svg", style=style, node_labels={}, mutation_labels={}, time_scale="rank"
)

[Schraiber]

This is super helpful! I think what's going on is that these are substitutions, since they seem to be occurring right at the type of the gene tree, as you can see in the attached figure. I think for my purposes I can just ignore these, since I'm ultimately interested in segregating mutations.