Plotting mutation rates over time

[hyanwong]

I would like to look at mutation rates over time, e.g. rates at different time slices, perhaps broken down by the different mutation types (transversions vs transitions, or triplet contexts).

I think the way to do this is to assign a mutational density to each edge (number of mutations on the edge divided by edge area), and then average the density of all edges that intersect a give timeslice, weighted by the proportion of the edge that exists in the slice.

This is a bit like the coalescence_rate calculation, I think. Has anyone done this before, and is there code to do so?

[nspope]

There's some code to efficiently get expected and observed mutational density in all inter-node intervals in tsdate (lacking an API: https://github.com/tskit-dev/tsdate/blob/3c08becf500dac3a422376126e70314064c360fe/tsdate/rescaling.py#L242). But, that's assuming mutations are all equivalent.

[hyanwong]

Neat, thanks! This does give me a start.

[hyanwong]

Here's what I have cooked up for the time being:

def mutation_rate(ts, timepoints):
    # return average mutation rates between timepoints
    mut_density = np.zeros(ts.num_edges)
    node_times = ts.nodes_time
    parent_times = node_times[ts.edges_parent]
    child_times = node_times[ts.edges_child]
    
    lengths = parent_times - child_times
    widths = ts.edges_right - ts.edges_left
    
    for m in ts.mutations():
        mut_density[m.edge] += 1
    mut_density /= (lengths * widths)
    
    # select regions of edges that intersect with a timeslice
    av_rate, mean_time = [], []
    for lo, hi in zip(timepoints[:-1], timepoints[1:]):
        lower = np.max([child_times], initial=lo, axis=0)
        upper = np.min([parent_times], initial=hi, axis=0)
        use = lower < upper
        weights = (upper[use] - lower[use]) * widths[use]
        av_rate.append(np.average(mut_density[use], weights=weights))
        mean_time.append((hi + lo)/2)
    return np.array(mean_time), np.array(av_rate)

Test:

ts = msprime.sim_ancestry(100, population_size=1e4, sequence_length=40e6, recombination_rate=1e-8)
ts = msprime.sim_mutations(ts, rate=1e-8)


plt.scatter(*mutation_rate(ts, np.append([0], np.logspace(2, np.log10(ts.max_time), 20))))
plt.xscale("log")
plt.title("vanilla msprime, mu=1e-8")
plt.xlabel(f"Time ({ts.time_units}")
plt.ylabel("Av mutation rate")

I assume the last few points are variable because they are based on very few edges / mutations