Repository navigation
Matching recombination rates for recapitation of SLiM with msprime #2333
Description
Hello, I am currently trying to simulate neutral and positive selection on a complex demography with multiple chromosomes.
My understanding is that for msprime, the recombination rate is the mean probability of a split event using a poisson distribution, so to simulate breakpoints for multiple chromosomes, the recombination rate at the breakpoint should be given as log(2).
For SLiM, my understanding is that it is the exact binomial probability, so the value should be given as 0.5.
I want to give a recombination rate in msprime 1e-08 within a chromosome, and this is where I am confused.
I wonder for which simulation should I give the value 1e-08 and for which simulation should I modify the value to match 1e-08 so that the two simulations agree in terms of recombination. My guess is that 1e-08 is the expected number of recombinations, so the msprime values should be 1e-08 and SLiM should be 1 - exp(-1e-08), which is anyway almost 1e-08 where the probability and the rate values are not so different in the first place, so it isn't much of an issue, but I just wanted to be sure.
Thank you for your time.
Activity
petrelharp commentedon Jan 8, 2025
Good question. I'll transfer this to a Discussion.
First, for between chromosomes, as you say, for msprime the value should be
log(2), and for SLiM the value should be 1/2 (see bits of the SLiM manual below).So, I think your question is for along chromosomes, when rates are very small (e.g.,
1e-8). Both software uses a Poisson model of recombination; the difference is that SLiM takes the user-provided per-bp probabilitypand transforms it to a rate byλ = - log(1 - p). So, if msprime has a given value ofλ, what value ofpshould we use? As you say,1 - exp(-λ). And, there is virtually no difference atλ=1e-8, since|exp(x) - 1 - x|isO(x^2):There are certainly bigger things to worry about in the simulation (e.g., discrete vs continuous time, bigger but still not worth worrying about).
Two caveats. First, note that since msprime's default model uses continuous time, we don't need to worry about "what if two breakpoints fall between two adjacent positions at once" like SLiM does. Second, I don't know why the
dtwfmodel doesn't need to worry about this (i.e., why islog(2)still the right answer; see the sorts of considerations in the SLiM manual). Probably whoever wrote that bit of the docs was considering this (@apragsdale?).