ms format sim_mutations() output

[mshpak76]

Is there a straightforward way to display the results of sim_mutations() in ms output format, i.e. a matrix with each row representing the genotype at every segregating site

For the msprime-generated object mts:

ts = msprime.sim_ancestry(samples={"A": 2, "B": 2, "C": 2}, sequence_length = 100000, demography = demography, recombination_rate = 1.14e-8, gene_conversion_rate=6.25e-8, gene_conversion_tract_length=1, random_seed=2222)
mts = msprime.sim_mutations(ts, rate = 1.14e-8, model = msprime.BinaryMutationModel(), random_seed = 1111)

I would like to display something along the lines of e.g. for 5 segregating sites

01100
01110
01001
01000
11100
11100

Presumably mts has the appropriate attribute to do so, but I don't know which one. Among the attributes listed is mts.genotype_matrix and mts.haplotypes, but I don't know if these are the correct attributes, nor what to do with them.

I would like to do this without running ms code in mspms.

Alternatively, my main reason for wanting this is to generate an output matrix from which I can readily compute allele frequencies at each segregating site and other population statistics, so perhaps there is a more efficient way of doing this.

[grahamgower]

Hey @mshpak76. The genotype_matrix() method of the tree sequence is probably what you want. There is also an (undocumented) tskit.write_ts() method that will give you ms-formatted output.

$ python -c "import tskit; help(tskit.write_ms)"
Help on function write_ms in module tskit.trees:

write_ms(tree_sequence, output, print_trees=False, precision=4, num_replicates=1, write_header=True)
    Write ``ms`` formatted output from the genotypes of a tree sequence
    or an iterator over tree sequences. Usage:

    .. code-block:: python

        import tskit as ts

        tree_sequence = msprime.simulate(
            sample_size=sample_size,
            Ne=Ne,
            length=length,
            mutation_rate=mutation_rate,
            recombination_rate=recombination_rate,
            random_seed=random_seed,
            num_replicates=num_replicates,
        )
        with open("output.ms", "w") as ms_file:
            ts.write_ms(tree_sequence, ms_file)

    :param ts tree_sequence: The tree sequence (or iterator over tree sequences) to
        write to ms file
    :param io.IOBase output: The file-like object to write the ms-style output
    :param bool print_trees: Boolean parameter to write out newick format trees
        to output [optional]
    :param int precision: Numerical precision with which to write the ms
        output [optional]
    :param bool write_header: Boolean parameter to write out the header. [optional]
    :param int num_replicates: Number of replicates simulated [required if
        num_replicates used in simulation]

    The first line of this ms-style output file written has two arguments which
    are sample size and number of replicates. The second line has a 0 as a substitute
    for the random seed.

[grahamgower]

It also occurs to me that because you're using a model converted from an ms command (as indicated by your previous discussion tskit-dev/msprime#2069), then your demographic model will correspond to haploid individuals. As such, you should specify ploidy=1 to sim_ancestry(). https://tskit.dev/msprime/docs/stable/ancestry.html#sec-ancestry-ploidy

This is all a bit confusing (I'm certainly confused), so we probably need a short vignette about it to help steer people in the right direction.

[grahamgower]

And to prove how confused I am, I just gave you wrong advice (the from_ms() docs in demes-python are wrong, we actually output sizes for diploid individuals). To match ms, you'll sample haploid genomes, but the demographic model is diploid.

ms_cmd = ... # something here
msprime.sim_ancestry(
    demography=msprime.Demography.from_demes(demes.from_ms(ms_cmd, N0=10000)),
    samples=[
        msprime.SampleSet(nsam, population=pop, ploidy=1)
        for pop, nsam in {"A": 2, "B": 2, "C": 2}.items()
    ],
    ploidy=2,  # this is the default
    ... # any other params
)

[jeromekelleher]

Alternatively, my main reason for wanting this is to generate an output matrix from which I can readily compute allele frequencies at each segregating site and other population statistics, so perhaps there is a more efficient way of doing this

It's probably worth while looking through a few of our tutorials @mshpak76, starting here. Tskit has lots of methods for computing popgen stats that will be much more efficient than outputting text.

[mshpak76]

I've been looking some of these over, and while there is a function to return the allele frequency spectrum, I'm unable to find a function that will (for k samples taken from k subpopulations) return the allele frequency at each segregating site for each of the k subpopulations. Is there some way to do this other than returning a matrix of haplotypes?

[jeromekelleher]

I don't think we have a function for doing that directly, but see #504 for some discussion and sample code.

To get the samples_sets for each of the populations, do samples_sets = [ts.samples(population=k) for k in range(ts.num_populations)]