Simulate diploid individuals using msprime

Question

[shahar]

Hi all,

I am using msprime version 1.1.1 with python 3.7.

I am simulating a relative simple simulation. The code:

POPULATION_SIZE = 1000
NUMBER_OF_SUBPOPS = 2
INDV_PER_POP = POPULATION_SIZE // NUMBER_OF_SUBPOPS

def run_simulation():
    pop_configs = [msprime.PopulationConfiguration(sample_size=INDV_PER_POP) for _ in range(NUMBER_OF_SUBPOPS)]

    ts = msprime.simulate(population_configurations=pop_configs, length=5e8, Ne=2000, mutation_rate=5e-7,
                          recombination_rate=1e-8,
                          demographic_events=[
                              msprime.MassMigration(10000, source=i, dest=0, proportion=1) for i in range(1, NUMBER_OF_SUBPOPS)
                          ])
    return ts

Later I want to output the simulation data as a VCF, using:

    with open(<some path>.vcf, 'w+') as f:
        ts.write_vcf(f)

I am surprised to see later that I have in my outputed VCF 1000 genomes. I asked for a sample size of 1000, and expected for diploid individuals. How can I simulate the same simulation with 1000 diploid individuals?

I was looking at the doc, and couldn't find where can I set ploidity,

Thanks a lot,
Shahar

Answer 1

[Ben]

You're using the legacy simulate function. It is highly recommended to use the updated methods of sim_ancestry and sim_mutations. There is a ploidy argument to sim_ancestry(). Have you also seen the section on ploidy in the docs which explains the effect of this in msprime? https://tskit.dev/msprime/docs/stable/ancestry.html#sec-ancestry-ploidy

[shahar]

Thanks, I think it is exactly what I looked for.