tree_sequence.genotype_matrix() contains non-segregating sites

[jeffspence]

Hello!

I would have thought that ts.genotype_matrix() would only contain segregating sites. But it seems like it sometimes contains monomorphic sites. Here's an example that generates this behavior:

import msprime as msp
import numpy as np


NUM_LOCI = 100_000_000
NUM_SAMPLES = 200
RECO_RATE = 1.25e-8
MUT_RATE = 1.25e-8

demography = msp.Demography()
demography.add_population(name='A', initial_size=1_000)
demography.add_population(name='B', initial_size=1_000)
demography.add_population(name='C', initial_size=1_000)
demography.add_population(name='D', initial_size=1_000)
demography.set_symmetric_migration_rate(['A', 'B', 'C', 'D'], rate=1e-3)

sampleset = [msp.SampleSet(NUM_SAMPLES // 4, population='A', time=0, ploidy=2),
             msp.SampleSet(NUM_SAMPLES // 4, population='B', time=0, ploidy=2),
             msp.SampleSet(NUM_SAMPLES // 4, population='C', time=0, ploidy=2),
             msp.SampleSet(NUM_SAMPLES // 4, population='D', time=0, ploidy=2)]

ts = msp.sim_ancestry(sampleset,
                      sequence_length=NUM_LOCI,
                      recombination_rate=RECO_RATE,
                      demography=demography,
                      random_seed=42)
mts = msp.sim_mutations(ts, rate=MUT_RATE, random_seed=42)
gts = mts.genotype_matrix()

print(gts.shape)
assert np.all(gts.sum(axis=1) > 0)

running the above on my MacBook gives that gts has shape (136359, 400) and then results in an assertion error.

I get this result running msprime==1.0.1 or msprime==1.0.0 on either my MacBook or on a linux server.

Another (minor) issue is that on my MacBook, gts has shape (136359, 400) but on the linux server it has shape (136808, 400) which seems weird since I believe I've seeded everything.

[benjeffery]

Hi @jeffspence thanks for posting your question.

It is expected that the genotypes contain non-segregating sites, as the simulation can add both back-mutations and silent mutations (as mentioned at https://tskit.dev/msprime/docs/latest/api.html#msprime.sim_mutations). The number of non-segregating sites should be consistent with the model and parameters.

As for differing results on differing OS installations, this is something that has occurred recently - we think due to differing versions of the GSL library that provides the randomness. How did you install msprime on each machine, pip, conda, or local build?

[jeffspence]

Thanks for this.

For genotype_matrix() -- it then includes the genotypes at all positions where there is one or more mutation even if back-mutations result in all individuals having the same genotype?

For GSL: both were installed using pip. One uses GSL version 2.6 and one uses GSL version 2.1

[petrelharp]

Said another way - genotype_matrix() gives you something that has (by default) dimensions equal to number of samples x number of sites, where the sites are positions at which there's at least one mutation. If you use msprime with the infinite alleles model and take the full genotype matrix then you won't get any monomorphic sites, but otherwise that's not guaranteed.

[petrelharp]

Perhaps there should be an option to return only segregating sites?

[jeffspence]

@petrelharp that would be nice! It would also be nice to get out just diallelic sites, since that's a somewhat standard pre-processing step for real pop-gen data (i.e., throw away sites with anything other than 2 alleles). Those are both fairly easy to do post hoc now, though, so it's no issue to me. I was just surprised by the non-segregating sites and wanted to make sure that that wasn't some sort of bug.

[andrewkern]

I agree @petrelharp, it would be nice to have an option to return only segregating sites

[petrelharp]

I agree that this is something we'd like to make easy for people. Maybe just a diallelic_only=False option to genotype_matrix( )? Hm, one problem with that is then that the rows of the genotype matrix are no longer in 1-to-1 correspondence with the sites, and so to connect the genotype matrix with any other information (e.g., matching the rows up with their positions along the genome) someone would then have to do some extra work.

So - is it really a common use case to want only the genotype matrix without the associated positions?

We have also floated the idea of adding a method that removes all monomorphic sites (although, simplify() does this already). So, I guess one solution is to simplify before outputting the genotype matrix. That wouldn't remove triallelic sites, though.

[hyanwong]

You could use delete_sites in conjunction with a helper function that identifies non-bialleleic sites, to make a new TS. The question is whether we want to provide such a helper function? If so, it could have options to only identify monomorphic sites, etc etc.

[petrelharp]

oh, great idea. It's on the list to write an alllele_frequency method; one might also ask for an array of the number of alleles, and that'd do it.