Basic question on recombination_rate and multiple chromosomes in msprime 1.0

[Richard Kerr]

Hi,

I have come back to revisit msprime after a long break and have discovered what I was used to is now Legacy (version 0.x) APIs.

In version 0.x I would use this (described by Jerome as a bit of a hack) to simulate multiple chromosomes (in this case 2 chromosomes)

rho=2e-8
positions=[0, 1.0e7-1, 1.0e7, 2.0e7-1]
rates=[rho, 0.5, rho, 0]
num_bp=int(positions[-1])

recombination_map = msprime.RecombinationMap(
    positions=positions, rates=rates, num_loci=num_bp)

I am not sure if version 1.0 has made it any less of a hack, but I am worried I am not doing it correctly. I gather one would now use 

rate_map=msprime.RateMap(position=positions, rate=rates)

and the positions and rates vectors are constructed a little differently. I now have
positions = [0.0, 9999999.0, 10000000.0, 19999999.0]

rates = [2e-08,  0.5,  2e-08]

this yields

┌────────────────────────────────────────────────────┐
│left      │right     │         mid│     span│   rate│
├────────────────────────────────────────────────────┤
│0         │9999999   │   4999999.5│  9999999│  2e-08│
│9999999   │10000000  │   9999999.5│        1│    0.5│
│10000000  │19999999  │  14999999.5│  9999999│  2e-08│
└────────────────────────────────────────────────────┘

I am not sure if the rates specified for any one span is the rate per base per generation within the span, or should I be specifying the rate over the span? Should rates for the first and last spans be 0.2 perhaps?

Hopefully someone can clarify this for me

RIchard

[Peter Ralph]

Hi, Richard! Thanks to discrete-sites recombination it's a lot less of
a hack, but you should probably read the section we have in the docs
on this:
https://tskit.dev/msprime/docs/latest/ancestry.html#multiple-chromosomes
especially:
https://tskit.dev/msprime/docs/latest/ancestry.html#example-2-simulating-without-a-pedigree

I think it will answer your questions (and post something if not!).

-- peter

> --
> You received this message because you are subscribed to the Google Groups "msprime-users" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to msprime-user...@googlegroups.com.
> To view this discussion on the web, visit https://groups.google.com/d/msgid/msprime-users/7148b2e1-a487-4b28-8e3e-263259f7baa2n%40googlegroups.com.

[Richard Kerr]

Thanks Peter,
that section certainly explains it. 

I am still mystified by the following logic:

"In the Poisson model of recombination used by msprime, the probability of co-inheritance across  generations of two adjacent base pairs separated by a recombination rate r per base pair is exp(-rt). To make these agree, we set the recombination rate in the base pair segments separating chromosomes to log(2.0)"

If there is any more literature that is able to explain this logic I certainly would benefit from it.

Kind regards,

Richard

[Peter Ralph]

Hm, which part do you want more explanation on? The "Poisson model"
part or why 2^{-t} is the right number, or what?

And - this seems like a good question that would be better answered as
a Discussion - mind posting it there?
https://github.com/tskit-dev/msprime/discussions

> To view this discussion on the web, visit https://groups.google.com/d/msgid/msprime-users/1d06d641-45c7-45ab-b4ea-f355cdfe903fn%40googlegroups.com.