Tracking IBD(Identity-by-descent) boundaries with and without full ARG

[jikhashkya]

Hi, I'm a new user to msprime and was playing around with how the tree sequences differed with and without record_full_arg option, specifically w.r.t to IBD boundaries. Since the number of trees reduce when disabling the option, I'm afraid this affects the intervals of the trees which in turn affects the IBD boundaries found for some sample pairs. I was advised to simplify() the ARG recorded tree sequence to see if IBD boundaries would be consistent, but I found that, that was not the case.
I used the following code for simulation:

num_samples = 4
ts = msprime.sim_ancestry(num_samples, sequence_length=10,
                          recombination_rate=0.1, random_seed=3,
                          record_full_arg=True, ploidy=2)
samples = ts.samples()
ts_simplified = ts.simplify(samples=samples)
ts_simplified.dump('data_dump-simplified.ts')
with open('pic-simplified.txt', 'w') as out:
    out.write(ts_simplified.draw_text())

I obtain the following tree for the code above:

Now, when I run the simulation with the same parameters but without the --record_full_arg option, I obtain the following tree:

This created inconsistency for IBDs found given the same simulation parameters. Here's an example for the IBDs found for sample pair (0, 4) for both trees:

• IBDs in the ARG recorded -> simplified tree sequence

SampleID-1	SampleID-2	Start Pos	End Pos	MRCA
0	4	0	3	18
0	4	3	5	16
0	4	5	6	15
0	4	6	7	17
0	4	7	10	15

• IBDs in the tree sequence without recording ARG

SampleID-1	SampleID-2	Start Pos	End Pos	MRCA
0	4	0	2	18
0	4	2	5	16
0	4	5	10	19

I'd love to know if this still makes sense and how it makes sense or, if that's not the case, how can the two cases be made consistent w.r.t IBD boundaries.

[petrelharp]

Hi, @jikhashkya! The problem here is that running the same code with and without record_full_arg=True does not get you the same tree sequence, even with the same random seed. In your example, you can see that the trees are in fact very simliar, but not equvalent. (By "equivalent" I mean having the same marginal genealogy at each position. On the other hand, these two things are the same:

import msprime, tskit
ts = msprime.sim_ancestry(4, sequence_length=10,
                          recombination_rate=0.1, random_seed=3,
                          record_full_arg=True, ploidy=2)
ts_simplified = ts.simplify()

a = ts.ibd_segments(within=[0, 4])
b = ts_simplified.ibd_segments(within=[0, 4])

Now, a and b are equivalent (happily):

>>> a
IbdResult({(0, 4): IbdSegmentList([IbdSegment(left=7.0, right=10.0, node=24), IbdSegment(left=5.0, right=6.0, node=24), IbdSegment(left=3.0, right=5.0, node=27), IbdSegment(left=6.0, right=7.0, node=31), IbdSegment(left=0.0, right=3.0, node=32)])})
>>> b
IbdResult({(0, 4): IbdSegmentList([IbdSegment(left=7.0, right=10.0, node=15), IbdSegment(left=5.0, right=6.0, node=15), IbdSegment(left=3.0, right=5.0, node=16), IbdSegment(left=6.0, right=7.0, node=17), IbdSegment(left=0.0, right=3.0, node=18)])})

Note that the nodes have different labels (since one tree sequence has more nodes than the other), but their times and IBD boundaries are the same.

Let me know if that's clear?

[petrelharp]

Right, that's what I was trying to say also. Why do you expect those two cases to be the same?

[jikhashkya]

I guess my assumption was that, just recording the ARG would not affect the IBDs existing in a simulated population if all the conditions and parameters were kept constant. Since both simulations (in theory) are the same , shouldn't the IBDs shared among individuals in both simulated instances remain the same irrespective of the detail of information stored in the tree sequences ?

[jeromekelleher]

I'm not clear on one detail here @jikhashkya - would you expect the IBD segments to be identical when running the two simulations (with and without record_full_arg) on the same seed, or do you expect the IBD segments to be identical in distribution (over many random replicates)?

As @petrelharp shows above, running a simulation with record_full_arg=True changes the random trajectory, so the output simulations will not be identical.

[jikhashkya]

Yes @jeromekelleher , I was expecting the IBD segments to be identical with the two simulations (with record_full_arg -> simplify() and without record_full_arg ) since in essence they are simulating the exact same population with the only exception being that I'm recording more details in the record_full_arg case. But yeah, my assumption didn't take into account the implementation specific details I suppose. Thank you all for the quick responses though, appreciate it.

[jeromekelleher]

It's pretty rare that you can depend on getting precisely the same output from msprime when varying parameters - the subtlest things can set the random trajectory off in a different direction. I've fallen into this trap many times!