Recording mutation frequencies through time

[leanngward]

I am working with simplified trees generated by SLiM. I have also overlaid neutral mutations using msprime.SLiMMutationModel.

I have been able to use the Python API to find the frequencies of mutations at the final generation (i.e., the time at which the tree was recorded):

sites = ts.sites()
for site in sites:
       tree = ts.at(site.position)
       for mutation in site.mutations:
               num_samples = tree.num_samples(mutation.node)
               total_samples = tree.num_samples()
               frequency = num_samples / total_samples

Question: Is there a way to find mutation frequencies at every generation in tree?

From what I could find in similar posts, I wonder if I might need to avoid simplifying trees in SLiM. However, I also noticed that the above code no longer returns correct frequencies when trees are not simplified, so I will have to make some adjustments. (Also, to clarify, "correct" frequencies are validated from mutation frequencies tracked during the SLiM simulation for non-neutral mutations.)

Thanks!

[petrelharp]

Question: Is there a way to find mutation frequencies at every generation in tree?

Depends what you mean here. If you mean "Is there a way to find mutation frequencies at every generation?", the answer is "no", because you've simplified, as you suspect. The way I'd do this is not by avoiding simplifying, but by Remembering all individuals alive at the times you want the frequencies at.

[leanngward]

That answer was helpful, thank you! I will share the code I used for tracking frequencies after making the changes you suggested (although it's probably not the most efficient way).

I did have a second question that arose from this one (and this may warrant a new discussion post, which I can do separately if needed).

I noticed that calculating frequencies as shown below no longer works when the simulated individuals have overlapping generations. It seems that the number of samples in the tree at each generation is not straightforward if individuals are older than one generation.

I found some resources for dealing with non-WF models on the tskit dev pages, but I did not see anything specifically addressing this. My follow-up question: Is it still possible to do mutation frequency calculation (per generation) for trees with overlapping generations?

Ncount = {}
NodeTimes = {}
MutCountPerGen = {}
MutFreqPerGen = {}
mutTypeDict = {}

nodes = tsnew.nodes()
for n in nodes:
        NodeTimes[n.id] = n.time
        if (n.time) in Ncount:
                Ncount[n.time] += 1
        else:
                Ncount[n.time] = 1


sites = tsnew.sites()
for site in sites:
tree = tsnew.at(site.position)
for mutation in site.mutations:
		time = math.floor(mutation.time)
		metadata = np.array(mutation.metadata["mutation_list"])
		muttype = metadata[0]["mutation_type"]
		mutTypeDict[mutation.id] = muttype
		while (time >= 0.0):
			if time in Ncount:
				mutcount = 0
				samplesundermutation = tree.samples(mutation.node)
				for m in samplesundermutation:
					if NodeTimes[m] == time:
						mutcount += 1
				key = str(mutation.id) + "," + str(time)
				MutCountPerGen[key] = mutcount
				MutFreqPerGen[key] = mutcount / Ncount[time]
				time -= 1
			else:
				time -= 1

[petrelharp]

I've not read your code in detail, but from what you say it sounds like a job for pyslim.individuals_alive_at(ts, ...) (also see some of the other bits of the pyslim documentation if you haven't already).

.... and okay, now I've read your code. I don't know what we get to assume about SLiM having put in the tree sequence enough to understand it. Have you seen the code in this thread? For instance, if there was only one mutation per site (or you were interested in only the "total derived allele frequency" per site) you could use this modified function (I think! not checked!!), for instance:

def allele_frequencies(ts, times):
    # times should be times at which SLiM has Remembered everyone
    sample_sets = [[n for ind in pyslim.individuals_alive_at(ts, t) for n in ts.individual(ind).nodes] for t in times]
    n = np.array([len(x) for x in sample_sets])
    def f(x):
       return x / n
    return ts.sample_count_stat(sample_sets, f, len(sample_sets), windows='sites', polarised=True, mode='site', strict=False, span_normalise=False)

It sounds like you're already doing some independent checking that you're getting the right answer. That's good, you'll need to do that.