Skip to content

Variable recomb - #11

Merged
jeromekelleher merged 4 commits into
tskit-dev:mainfrom
jeromekelleher:variable-recomb
Oct 28, 2026
Merged

jeromekelleher merged 4 commits into
tskit-dev:mainfrom
jeromekelleher:variable-recomb

Conversation

@jeromekelleher

Copy link
Copy Markdown
Member

Port variable recombination infrastucture from tsinfer

Copy read_genetic_map and recombination_probabilities from tsinfer
(tsinfer.matching, commits 01e3aed and 16c4567), which convert a constant
rate per bp or a HapMap genetic map into per-site recombination
probabilities with Haldane's map function, with site 0 = 0. These are the
r_j that tsls.Hmm takes: tsls interprets recombination_rate[j] the same
way as tsinfer (recombining before site j), so the probabilities can be
passed straight through.

The tests are tsinfer's, with msprime's RateMap as an oracle, plus three
tsls.Hmm checks: a high rate gives a switch, a low rate with a high
mismatch probability gives a mismatch instead, and a map with zero rate
between two sites forces the switch elsewhere. msprime is added to the
test dependency group.

tsinfer's HmmParameters and config plumbing are not ported: they are
specific to tsinfer's pipeline.
run.py takes --recombination-rate (a rate per bp, or the path to a
HapMap genetic map) and --mismatch, which replace the model.json rho and
mu as tsinfer's recombination_rate and mismatch would: the per-site
recombination probabilities come from tsls.model.recombination_probabilities
at the reference's sites. The same workload, whose truth paths tsinfer
found with constant rho = 0.01, can then be rerun under a map.

- tsinfer-hmm and tsls engines take rho as a scalar or a per-site array.
- Results have a recombination_rate column, and rho is the mean over
  sites 1.. when it is per-site. run.py refuses to append to a results
  file with the old columns, rather than misaligning rows.
- tsinfer_data.load_model reads recombination_rate from a tsinfer config
  (rejecting it with recombination, as tsinfer does), and run.py uses it
  from model.json.
- tsinfer_summary.py reports the recombination and mu of each run.
genetic_map.py writes a HomSap genetic map from stdpopsim (by default
HapMapII_GRCh38 chr20) as a HapMap file, and checks that
tsls.model.read_genetic_map reads back msprime's cumulative mass at
100k positions. stdpopsim is not a dependency; run with
uv run --with stdpopsim.

The chr20 map has 82,979 intervals and 1.081 Morgans over 64.4 Mb. On the
580,909 chr20p sample inference sites (0.23-26.3 Mb, 0.508 Morgans), the
per-site probabilities are identical to those from tsinfer's
matching.recombination_probabilities. Their median is 8.4e-8 (10% 2.8e-9,
99% 1.6e-5, max 6.2e-4, none zero), against tsinfer's constant 0.01: with
n = 432,853 nodes the median switch odds are 1.9e-13 against 2.3e-8.
Rerun the chr20p samples workload (200 queries) with per-site
recombination probabilities from the HapMap II GRCh38 chr20 map in place
of tsinfer's constant rho = 0.01, run 4 at a time on 4 cores:

  recomb  mu     engine   queries segs  mm    time  MiB  traceback
  0.01    1e-20  tsls     200     244   0     2.2s  12   13
  map     1e-20  tsinfer  33      159   20.7  12.2s 349  113
  map     1e-20  tsls     30      166   21.6  4.5s  36   106
  map     1e-8   tsls     200     130   86.8  4.6s  37   114
  map     1e-4   tsls     200     108   130   4.9s  38   134

No query keeps tsinfer's path and only 2-3% of its breakpoints survive.
Mismatches appear even with mu = 1e-20: under the map two switches (away
and back) have odds of about 4e-26, less than one mismatch. Matching is
2-5x slower with 3x the memory and 8x the traceback, because smaller
switch odds mean fewer equal likelihoods to merge; tsinfer-hmm's memory
grows far more than tsls's. tsinfer-hmm and tsls give identical paths on
the 30 queries both ran.

The map run with mu = 1e-20 was stopped early (tsinfer-hmm at about
12 s per query); the results CSVs are not tracked.
@jeromekelleher
jeromekelleher merged commit 3447ab7 into tskit-dev:main Oct 28, 2026
28 of 30 checks passed
@jeromekelleher
jeromekelleher deleted the variable-recomb branch October 28, 2026 16:36
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant