using psmc on msprime output

[alanrogers]

To run msprime output through psmc, I need to produce that output either in fastq format or in psmc's psmcfa format. Are there examples of how to do this?

Thanks, Alan

[jeromekelleher]

No, I'm afraid not. The simplest approach is probably to output to vcf or FASTA and convert from there?

[grahamgower]

There's some code to do this in the xsmc repository that should be very easy to pull out and use. https://github.com/terhorst/xsmc/blob/e0f0e558590f4b10bd1eef45567c10afc34dde43/xsmc/supporting/mspsmc.py#L38-L54

[grahamgower]

Oh, and there's https://github.com/terhorst/mspsmc

[alanrogers]

The mspsmc repo looks ideal, but the documentation is a broken link.

[alanrogers]

I'm experimenting with the code in the xsmc repo.

[alanrogers]

For the benefit of anyone else trying to do this, here is what I came up with, based on the xsmc code. chr_generator is the generator returned by msprime.sim_ancestry. I was using the num_replicates argument. If you don't use that argument, the return type is different. chr is the current replicate (representing a chromosome) and i is its index. u_per_site is the mutation rate per base pair.

    # Put mutations onto the gene genealogy and output in psmcfa format
    w = 100
    for i, chr in enumerate(chr_generator):
        sim = msp.sim_mutations(chr, rate = u_per_site)
        L = int(sim.get_sequence_length()) // w
        outstr = ["T"] * L
        for v in sim.variants():
            gt = v.genotypes
            if gt[0] != gt[1]:
                outstr[int(v.position) // w] = "K"
        print("> chr%d" % i)
        print("\n".join(textwrap.wrap("".join(outstr), width=79)))
        print()

Thanks for the help.