Edge_Diffs() with >4 edges_out (and edges_in)

[mebh]

Hi, I'm big fan of the tskit and msprime software, use them all the time. In the course of a recombination inference analysis, I noticed that for some tree transitions there are more than 4 edges_out (= #edges_in). As far as I can figure, a single recombination event should only be compatible with 2,3, or 4 edges_out (= #edges_in). I am a little confused as to how this makes sense. It only occurs for a small fraction of transitions for vanilla simulations of a single pop. In msprime, can a single transition possibly encode multiple recombination events? Any pointers as to what is possibly going on?

As far as I can tell this does not happen for a vanilla ms.simulate() tree_seq, but it does happen when using ms.sim_ancestry() for what I think should be the same single pop, constant size, simulation. See examples below to see what I mean exactly.

The # edges in/out = 2,3,4 for a single recomb event seems to accord with Kelleher, Jerome, Alison M. Etheridge, and Gilean McVean. "Efficient coalescent simulation and genealogical analysis for large sample sizes." PLoS computational biology 12.5 (2016), but I am not sure on that.

Thanks!

Example cases
ms.simulate:

ts = ms.simulate(sample_size=30, Ne=1e4, recombination_rate=1e-8, mutation_rate=1e-8,length=5e6, )
counts = Counter([len(edges_out) for (interval,edges_out,edges_in) in ts.edge_diffs()])
print(counts.most_common())
>[(4, 3817), (2, 1387), (3, 1161), (0, 1)]

ms.sim_ancestry:

demography = ms.Demography()
demography.add_population(name="Pop1", initial_size=1e4)
ts = ms.sim_ancestry(
        recombination_rate=1e-8,
        sequence_length=5e6,  
        samples=[
            ms.SampleSet(30, ploidy=1, population='Pop1'),
        ],
        demography = demography,
        random_seed=1001
    )
counts = Counter([len(edges_out) for (interval,edges_out,edges_in) in ts.edge_diffs()])
print(counts.most_common())
>[(4, 3903), (2, 1438), (3, 1310), (7, 2), (8, 2), (0, 1), (6, 1)]

[petrelharp]

Good question, @mebh! With msprime v1.0, we've switched to a default of a discrete genome, i.e., having recombinations occur at integer positions. This applies to sim_ancestry, as you've noted, but not to simulate, as we wanted to keep things backwards compatible. So, what you've observed is just what I'd expect - a small number of breakpoints that have been hit multiple times through the course of the simulation, just as now by default sim_mutations will produce a (usually small) proportion of sites with multiple mutations.

[mebh]

Got it, thanks @petrelharp! Yeah I wondered if it was something like that- I dimly recalled a discussion of the discrete v continuous in the (old?) docs but couldn't find it again. And I now see that the interval positions are indeed all at integer positions for the ms.sim_ancestry() trees while they are reals for ms.simulate(). So that all check out :) Thanks for clearing that up for me!