How to remove low frequency mutations?

[kliao12]

Hi,

I'm currently trying to subset the mutations in a simulated tree sequence based on allele frequency. When trying to run very similar code as the one provided in the tutorial (copied below), I'm getting the error "AttributeError: 'TreeSequence' object has no attribute 'set_mutations'". It seems the set_mutations function was removed at some point (if I install an older version i.e. 3.6.2 set_mutations exists). In the latest versions of msprime, is there an alternative to subset mutations? Apologies if I'm missing an obvious alternative function.

        import msprime
        
        def set_mutations_example():
            tree_sequence = msprime.simulate(
                sample_size=10000, Ne=1e4, length=1e7, recombination_rate=2e-8,
                mutation_rate=2e-8)
            print("Simulated ", tree_sequence.get_num_mutations(), "mutations")
            common_mutations = []
            for tree in tree_sequence.trees():
                for mutation in tree.mutations():
                    p = tree.get_num_leaves(mutation.node) / tree.get_sample_size()
                    if p >= 0.5:
                        common_mutations.append(mutation)
            tree_sequence.set_mutations(common_mutations)
            print("Reduced to ", tree_sequence.get_num_mutations(), "common mutations")

Thanks,
Kevin

[petrelharp]

Gee, thanks for the report! It's true, that method isn't around any more (it was never documented, so wasn't supposed to be something that we supported long-term). How did you find that tutorial, btw? I found it in this old version of the docs by googling, but I'm hoping others won't come upon this old version very easily.

If you want to do that now, here's what I'd do:

import msprime, tskit

tree_sequence = msprime.simulate(
    sample_size=10000, Ne=1e4, length=1e5, recombination_rate=2e-8,
    mutation_rate=2e-8)
print("Simulated ", tree_sequence.get_num_mutations(), "mutations")

tables = tree_sequence.tables
tables.mutations.clear()
for tree in tree_sequence.trees():
    for mutation in tree.mutations():
        p = tree.get_num_leaves(mutation.node) / tree.get_sample_size()
        if p >= 0.5:
            _ = tables.mutations.add_row(
                 site=mutation.site,
                 node=mutation.node,
                 derived_state=mutation.derived_state,
                 metadata=mutation.metadata,
                 time=mutation.time)

tables.compute_mutation_parents()
ts = tables.tree_sequence()

print("Reduced to ", ts.get_num_mutations(), "common mutations")

Note that this will still have invariant sites; if you'd like to remove these then you could run .simplify(), for instance.

If there's a bug here, it's that you found this old documentation for something that was never supposed to be released. But, it's a good question that probably belongs in Discussions - I'm going to transfer it over there (and let me know if you don't agree!). Thanks again!

[benjeffery]

Hi @kliao12!

(seems @petrelharp and I were typing at the same time!)

set_mutations was removed way back in 2017 in before tskit was split from msprime 805d745

The tutorial example was removed in 2018: 5b25db2 which would have been released in 0.5.0 in Feb '18

Where were you viewing the tutorial? It would be good to know so we can update it!

The TreeSequence class is now immutable, if you want to modify the mutations table you create a new one via the tskit.TableCollection which provides ways to mutate the underlying data of a tree sequence:

tables = tree_sequence.tables
tables.mutations.clear()
for tree in tree_sequence.trees():
    for mutation in tree.mutations():
        p = tree.get_num_leaves(mutation.node) / tree.get_sample_size()
        if p >= 0.5:
            tables.mutations.append(mutation)
tables.compute_mutation_parents()
ts = tables.tree_sequence()

(@petrelharp used add_row but append is the new way since 0.3.6 :D )

[petrelharp]

Whoops; couldn't remember if append had been released. =) This was is nicer!

[hyanwong]

I assume this worked for you @kliao12 ? I've marked it as the "answer"

[benjeffery]

Read the docs have also removed the old tutorial from a neglected cache.

[kliao12]

Great thanks to both of you for the prompt responses and suggestions on working with the now immutable TreeSequence class! And yeah I was working off the apparent old docs (https://msprime.readthedocs.io/en/multi-requirements/tutorial.html), which was the first link that popped up for me when googling "msprime tutorial".

[benjeffery]

Thanks for letting us know about that site, I thought I had disabled all the old docs. I've raised an issue to fix this at #1691

The tutorials for the tree sequence ecosystem are now at https://tskit.dev/tutorials/,