tskit-dev / tskit

Population-scale genomics
MIT License
155 stars 73 forks source link

`pair_coalescence_rate` doesn't account for empty regions (e.g. flanks) #3053

Open hyanwong opened 6 hours ago

hyanwong commented 6 hours ago

You might hope that the span-normalised coalescence rate would be the same for a given tree, even if the tree sequences has empty regions at the start or end, or even in the middle. However, this test fails:

ts = tskit.Tree.generate_comb(10, span=10).tree_sequence
time_intervals = np.logspace(0, np.log10(ts.max_time), 4)
time_intervals = np.concatenate(([0], time_intervals, [np.inf]))
ts_with_missing = ts.delete_intervals([[0,2],[6,8]])
assert np.allclose(
    ts.pair_coalescence_rates(time_intervals),
    ts_with_missing.pair_coalescence_rates(time_intervals),
)

It passes once we trim down:

ts_trimmed = ts.delete_intervals([[0,2],[6,10]]).trim()
assert np.allclose(
    ts.pair_coalescence_rates(time_intervals),
    ts_trimmed.pair_coalescence_rates(time_intervals),
)

@nspope says

I see what the problem is. Basically, I'm assuming weights sum to 1 so that we can calculate rates even when there are multiple roots (e.g. decapitated trees). That won't be the case if we're incorrectly normalising, however

nspope commented 6 hours ago

This seems like a pretty serious usability issue, as there will often be interior regions that are masked (e.g. geneic regions) which can't easily be trimmed off. This is an easy fix, I'll get it in ASAP.

hyanwong commented 5 hours ago

Thanks @nspope . Note that the other stats methods also have a similar (but not identical) issue: https://github.com/tskit-dev/tskit/issues/287