Msprime: [Docs] Update section on "simulating multiple chromosomes"

[jeromekelleher]

This is currently out of date and, should be using the discrete option to RecombinationMap.

Also some important details should be spelled out. @apragsdale, would you mind listing out what needs to be said here?

[mufernando]

I think I could help with this, but it's not clear to me what should be done. @apragsdale ?

[jeromekelleher]

Thanks @mufernando, that would be really helpful. Just to note also that the bit above is also out of date now, and we should be using the discrete_genome option to simulate.

[apragsdale]

Thanks @mufernando. I haven’t kept up with all the changes in the last few months that have been going on - so I need to see how the discrete recombination works as well... not sure, but it’s possible it now allows for point probabilities if recombination between loci?

I’ll get myself up to speed, hopefully soon, and would be able to help out with this.

[petrelharp]

We can now do recombination with discrete (i.e., integer) coordinates, so I think all that needs to be said is that you should put a recombination rate of 0.5/bp on a 1bp gap between chromosomes? And be careful to get the endpoints right?

[apragsdale]

That sounds right to me. What I was just trying to check, looking through old closed issues and PRs, is would you get a huge number of trees in between the adjacent loci where you have recombination rate of 0.5? (Or is it now a recombination probability, which is different?) I didn't see if or where that was explicitly described, though.

[petrelharp]

No - you should (only!) do this with discrete recombination, in which case recombination breakpoints only occur at integer positions.

[petrelharp]

And, note that discrete recombination will be the default for the new simulation interface. (but we should still say explicitly what's going on)

[apragsdale]

Ok great!! Thanks Peter. And in the discrete genome setting, a recombination rate of 0.5/bp between those loci really does mean "recombination event occurs with probability 0.5" instead of "number of recombination events are drawn from Poisson(0.5)" (and then you sum over probabilities of odd numbers of recombinations), right? If that's the case, then perfect!

Another way of simulating multiple chromosomes is running independent simulations through a pedigree. However, I don't know what the state of the pedigree-based simulations is at the moment. But the benefit in the pedigree setting is that the different chromosomes can be simulated in parallel, instead of in series in the same simulation with a long recombination map that has been spliced together. The independent simulations through the same pedigree imposes recombination probability of one half between all those parallel-simulated segments.

[petrelharp]

We're by default in continuous time, so it means that "recombinations occur at a rate of 0.5 per unit time" (and I guess we should note that if a unit of time is not a generation they should adjust accordingly). I think this is correct, because it gives the correct long-run density of recombination events across time. I don't know what the DTWF does, although I know that migration rates are actually migration probabilities, and so if recombination is similar, that is correct also. I guess someone should look at the dtwf code to confirm this.

[petrelharp]

Ah, good point about the pedigree! But I don't think the pedigree simulation stuff is really ready?

Huh - I've thought in the past that running different chromosomes independently is equivalent to running them at the same time with a 0.5 recombination bit between them, in the Hudson coalescent. But right now that's not seeming correct to me. Is it?

[apragsdale]

No, I don't think it is ready. But I'll check in with Dominic and Ivan who were the latest to be working on it, I think.

In the Hudson coalescent, it isn't quite correct, because the recombination events split those segments into more than two parental lineages. But that is correct in the DTWF model, because multiple recombination events along the genome bring you back to the same parental genome.

[petrelharp]

In the Hudson coalescent, it isn't quite correct, because the recombination events split those segments into more than two parental lineages. But that is correct in the DTWF model, because multiple recombination events along the genome bring you back to the same parental genome.

I'm confused about what "that is correct" is referring to, but, to clarify: it is definitely not true that simulating different chromsomes in completely separate runs of msprime under the DTWF is equivalent to simulating those chromosomes together in the same run, right? I was asking whether these two things were equivalent: (a) simulating chr1 and chr2 in separate runs of msprime(model='hudson'), then putting them together, and (b) simulating chr1+chr2 in a single run of msprime(model='hudson').

[apragsdale]

Sorry, I wasn't very clear. I think for simulating multiple chromosomes, Hudson is never really the right choice, if you care about getting the long-range correlations in the gene genealogies correct. Imagine simulating chr1+chr2+chr3+... with 0.5 between each. In the Hudson model, each recombination event takes you to a new parent, so you could end up with many more than two parents for a given genome. Then chr1 and chr3 will be less correlated than chr1 and chr2 or chr2 and chr3.

Simulating multiple chromosomes appended together with the DTWF model avoids this issue, because multiple recombination events just alternate between the two parents, so all pairs of chromosomes are equally correlated on average (whether they are neighbors or whether you are looking at chr1 vs chr22).

In the pedigree setting, you are enforcing a DTWF type model (without having to specify it) because there are always two parents per individual in the pedigree. Independent simulations through the same pedigree will result in the same desired correlation between chromosomes as if you simulated them all appended together under the DTWF model (but not the Hudson model).

[petrelharp]

Ok, great. This is a good summary of the issues. (And, the answer to my question is "no, those aren't equivalent", for these reasons.) I think we should not alarm people too much about Hudson, since it's pretty darn close. And since dtwf is substantially slower (I think this is still true?). But we could recommend a final period of dtwf - do you think that's a good idea? Do you know a good recommendation of how long this period should be?

[apragsdale]

Right, sorry, I didn't answer your question directly (though I meant to :-) ). I think you're right, that for most use cases Hudson is good enough - though in those cases, where you're not looking at long-range statistics like LD between distant loci or IBD between close relatives, for example, you probably don't gain anything by stitching chromosomes together. It just ends up adding a lot to run time (quadratic in genome length, right?), and doesn't meaningfully change the statistics on the tree sequence that you care about.

If you do care about correlations across chromosomes, or care about close relatives, you'd want DTWF but I don't think you need many generations of DTWF.

In Dominic's paper describing the DTWF (https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1008619), we highlight how using DTWF vs Hudson affects distribution of IBD for relatives (Fig 3) and variance of ancestry proportions after admixture (Fig 4). From Fig 4, you could guess that you'd do a good enough job with 10 generations of DTWF. Also in that paper, long-range LD (Appendix S2, see panel B) looks to need even fewer generations of DTWF to look right. So my non-rigorous recommendation would be 10 generations to be "close enough", and maybe 20 generations to be on the safe side. 10 or 20 generations also has a pretty negligible effect on the run time, I'd think.

[jeromekelleher]

Great discussion, thanks @apragsdale and @petrelharp!