removing fixed mutations

[petrelharp]

We can at times end up with a large number of fixed mutations. It would be nice to have a way to remove these. I think we discussed having simplify() do this, but it does not:

ts = msprime.simulate(10, recombination_rate=2, mutation_rate=2)
sts = ts.simplify([0,1])
sts.genotype_matrix()
### 
array([[1, 1],
       [1, 1],
       [1, 1],
       [1, 0],
       [1, 1],
       [1, 1],
       [1, 1],
       [1, 1]], dtype=uint8)

This would require a bit of bookkeeping (e.g. updating the ancestral state at sites that still have a segregating mutation). Not sure whether to propose this as an argument to simplify (filter_fixed_mutations) or as its own function. This is not urgent, as far as I know, but if we were to change the default behavior of simplify, it would be nice to do it soon.

[jeromekelleher]

I remember going through this, and it turns out to be quite tricky in the general case of lots of mutations at a site to know that it's fixed. I think the semantics that we've ended up with are good, in that we provide the option to remove objects that have no references (sites, populations, etc). Fixed sites are quite a different thing.

I think we can do it with a function easily enough though? Should be fast enough in Python, since we can do in in O(#mutations) per site, using the counting logic as we have in the site general stats algorithm.

[petrelharp]

Here's a python version, is this what you had in mind?

def remove_monomorphic(sts)
  tables = sts.tables
  tables.sites.clear()
  tables.mutations.clear()
  n = sts.num_samples
  for t in sts.trees(sample_counts=True):
     for s in t.sites():
        visited = False
        for m in s.mutations:
           k = t.num_samples(m.node)
           if k > 0 and k < n:
               if not visited:
                   visited = True
                   sid = tables.sites.add_row(s.position, s.ancestral_state, metadata=s.metadata)
               tables.mutations.add_row(sid, m.node, m.derived_state, parent=-1, metadata=None)
  tables.compute_mutation_parents()
  return tables.tree_sequence()

nts = remove_monomorphic(sts)
nts.genotype_matrix()

###
array([[1, 0],
       [0, 1],
       [1, 0],
       [0, 1]], dtype=uint8)

[jeromekelleher]

Basically, yeah. But this won't handle funky situations where we have (say) n mutations to 1 over the leaves, will it? But, we should be able to build an allele count table like we do in general_site_stats and reason using that, right?

[petrelharp]

True, but I don't think we want to remove those mutations.

[jeromekelleher]

Really? The site is still monomorphic right? What's the difference?

[petrelharp]

Well, in that case everyone is IBS but not IBD. They're different cases, even if you can't tell from the sequence.

[jeromekelleher]

I see that, but it's a statement about the sequences we're making, not the ancestry I would have thought. The semantics will get weird and strained if we don't I think (remove_monomorphic_if_single_mutation).

[petrelharp]

So, I don't mind if there's a few fixed mutations in there, especially if they are due to multiple mutations at the same site. The main reason I wanted this is because after a long SLiM simulatoin, a large portion of the tree sequence file can be due to fixed mutations, even after simplification (they are still there because of the initial generation). So, this is purely pragmatic, about file size. But, I don't want this to discard the still-segregating-but-invisibly multiple IBS mutations, because that is discarding information relevant to polymorphism in the population. I can imagine someone wanting the method that you're talking about, too, but I don't think it's important because a bit of bookkeeping can be used to ignore those.

[petrelharp]

For instance, with this SLiM recipe (100 individuals for 100,000 generations):

initialize()
{
    setSeed(23);
    initializeTreeSeq();
    initializeMutationRate(1e-8);
    initializeMutationType("m1", 0.5, "f", 0.0);
    initializeGenomicElementType("g1", m1, 1.0);
    initializeGenomicElement(g1, 0, 1e8-1);
    initializeRecombinationRate(1e-8);
}

1 early() {
    sim.addSubpop("p1", 100);
}

100000 {
    sim.treeSeqOutput("fixed.trees");
    catn("Done.");
    sim.simulationFinished();
}

we see that 97.6% of the mutations are fixed:

import pyslim

ts = pyslim.load("fixed.trees")
sts = ts.simplify()

freqs = ts.general_stat(np.ones(ts.num_samples).reshape((ts.num_samples, 1)), lambda x: x, windows="sites", polarised=True, span_normalise=False).reshape((ts.num_sites,))

np.mean(freqs == 200)

... and in Anastasia's simulations, this makes the metadata column longer than is allowed.

ps. we should make an "allele frequency" function, which we can do as above.

[jeromekelleher]

Ah, I see. Should we add a filter_fixed_simple_sites (or something) option to simplify then? I think catching the single mutation fixed sites is pretty easy? I don't want to call it monomorphic or whatever as I think this is misleading in the multiple mutation case.

[jeromekelleher]

ps. we should make an "allele frequency" function, which we can do as above.

+1

[petrelharp]

Just filter_fixed_mutations should do it? (i.e., fixed mutations, not fixed alleles...)

[jeromekelleher]

Just filter_fixed_mutations should do it? (i.e., fixed mutations, not fixed alleles...)

... Yeah, OK. But this implies that we only remove the mutation, not the site. However, if this is applied before filter_sites, we'd remove the site because any fixed mutations would have already been filtered out.

By fixed mutation, we mean there is exactly one root and there is a mutation over it. What if there are multiple mutations over the root? Do remove all of them, or do we leave them be? (Seems best to remove all of them?)

[petrelharp]

Hm, well we don't necessarily remove the site, since what if there's a polymorphic mutation at it still...

Do we remove all of them, or do we leave them be? (Seems best to remove all of them?)

All of them!