Stdpopsim: Drosophila two population msprime simulations slow with recombination map

[tjstruck]

I've been trying to add drosophila for the two population analysis, however it seems that incorporating the recombination map into the msprime simulations is slowing it down a lot (ex. I have a simulation with 2 samples per population running and has been running for >5 days).

Is this expected or is there an error on my part? See zipped folder with a standalone script (dmel_two_pop_recomb_example.py) and scripts needed to run snakemake (edited to just run msprime).

dmel_two_pop_recomb_test.zip

Let me know if there is any other information that may be useful.

Update:
It took 5 to 6 days for the msprime simulations with two samples per-population to finish.

[jradrion]

Andy and I have been looking into this, and we think this is standard behavior and not a bug. I had previously used msprime to simulate the same dmel recombination maps under equilibrium for a different project, and I was experiencing run times of ~1hr for simulating a single chromosome (2R) with Ne =2.5e5. The human models take significantly longer to run as well if you bump Ne up > 100,000. It might be worth considering abandoning recombination maps for models with very large Ne, if we will need to do many simulation reps.

[jeromekelleher]

Having a big Ne really pushes up the effective recombination rate, so it does make simulation times long.

[jradrion]

I've been running a few tests to check the run times for simulating only chr2R under Li & Stephan (population 0) for samples sizes 2,4,8,16,and 32. Sample sizes 8, 16, and 32 are still running, but since n=2 and n=4 have now finished running, I figured this would be worth sharing.

The run times were 4718 and 8096 minutes, for n=2 and n=4 respectively.

We should probably be thinking about exactly what populations we want to ultimately sample (Africa, Europe, both) and what those samples sizes will be, so that we don't have to simulate more than once. I can then set these up to run on our machine in the Kern lab, and we can just let it go for as long as it takes.

[jeromekelleher]

Jeez that's slow. Must be a massive recombination rate... I wonder how realistic this is, and whether the issues we're talking about in this preprint are relevant. Might be worth trying this with a few 1000 generations of DTWF to see what happens? Or is this opening up another can of worms?

[andrewkern]

what do you think @jeromekelleher - should this model be feasible and we have a bug somewhere?

[jeromekelleher]

what do you think @jeromekelleher - should this model be feasible and we have a bug somewhere?

Looks fine to me. I'm just wondering if all the recombination is actually meaningful - how many ancestors are we getting per generation in the recent past? Might be worth trying with DTWF just to see if there's much difference? Just involves setting model="dtwf" when calling simulate.

[andrewkern]

yeah we can give the DTWF a whirl for sure.

[jeromekelleher]

Worth a go - should speed things up in principle, and be more realistic.

[jradrion]

Sounds good. I'm on it!

[jradrion]

Just an update: the DTWF model for Li & Stephan (n=2; chr2R) is still running 14 days later.

[jradrion]

Runtimes for the other tests are as follows:

chrom | simulation_model | sample_size | demographic_model | time (min) |  
-- | -- | -- | -- | -- | --
2R | dtwf | 2 | LiStephan |  ? |  
2R | standard_coalescent | 2 | LiStephan | 4718 |  
2R | standard_coalescent | 4 | LiStephan | 8096 |  
2R | standard_coalescent | 8 | LiStephan | 10506 |  
2R | standard_coalescent | 16 | LiStephan | 12213 |  

[jeromekelleher]

Thanks @jradrion - looks like that's not a silver bullet then. :disappointed:

[jradrion]

Unfortunately, no.

[jradrion]

Update: Hybrid sims don't appear to be any faster in this instance, so I've moved ahead with the genome-wide standard coalescent sims. I'm currently running 10 reps, 50 samples, for both the African and European pops using LiStephanTwoPopulation. We have a walltime limit of 30 days; I'm crossing my fingers that they finish in time. I will provide the tree sequences once they have been generated.

[jradrion]

@everyone I have 10 reps of LiStephanTwoPopulation simulations for those that would like them.

Available here

These are using n = 50, but only simulating a chunk of chr2R (chr2R:10000000-13000000). The whole genome simulations are still running.

[ckyriazis]

Thanks @jradrion! I will get working on producing preliminary results for the analysis pipeline.

[jeromekelleher]

I think we can close this --- these sims are slow and that's just how it is I think.

[ckyriazis]

@jradrion any updates on the whole genome simulations?

[jradrion]

@ckyriazis Yeah, they were recently killed due to cluster maintenance after running for at least 30 days. I think we're going to need to re-scale the population sizes. I'm messaging Andy about this now. I'll get back to you ASAP.

[ckyriazis]

Ahhh bummer, let me know what he says.