read-in custom ratemap

[archgen]

Hi, just checking if it is possible to read in my own ratemap.txt file in the ratemap format of (left \t right \t mid \t span \t rate)? I see there is the read_hapmap() but wondering if msprime can read-in a custom ratemap.
Thanks very much!!

[grahamgower]

Hey @archgen! 👋

This should be possible, assuming your rate column uses the units cM/Mb (if not, you'll have to convert), and assuming your map intervals are contiguous (otherwise you'll have to figure out what rate to give the "gaps"). Something like the following should work.

chromosome_length = ... # some appropriate number here
ratemap = msprime.RateMap.read_hapmap(
    "ratemap.txt", position_col=0, rate_col=4, sequence_length=chromosome_length
)

If you get errors from this, you might need to manually edit your file (e.g. add entries to ensure the telomeres have rate zero). See the API docs for info on the meaning of the columns.
https://tskit.dev/msprime/docs/stable/api.html#msprime.RateMap

And then see the discussion of rate maps to see how to inspect the data structure once you've read it in, to double check that you got the ratemap you expect.
https://tskit.dev/msprime/docs/stable/rate_maps.html

[jeromekelleher]

Another option is to use numpy.loadtxt to get the columns of your file as numpy arrays and pass these to the RateMap constructor directly (this is what read_hapmap is doing under the hood, with some extra complications for dealing with quirks in the input)