Using deCode rate map not producing recombination events

[jikhashkya]

Hi,

I was trying to use the deCode rate map to introduce non-uniform recombination into a simulation using the following piece of code:

import msprime
import tskit
import numpy as np

num_samples = 4000
decode_map_file = '</path/to/decodemap>'
rate_map = msprime.RateMap.read_hapmap(decode_map_file, map_col=3)
ts = msprime.sim_ancestry(
    samples=num_samples,
    population_size=1,
    ploidy=2,
    random_seed=123,
    recombination_rate=rate_map,
    record_full_arg=True
)
ts = msprime.sim_mutations(ts, rate=0.001, model=msprime.BinaryMutationModel(), random_seed=1234)

re_nodes = np.where(tables.nodes.flags == msprime.NODE_IS_RE_EVENT)[0]
if re_nodes.size <= 0:
    print("No recombination nodes ???")
    exit(1)

And the deCode map's first few lines are as follows:

Chromosome      Position(bp)    Rate(cM/Mb)     Map(cM)
chr21   14180791        0.0     0.0
chr21   14181208        6.112705354778751e-25   2.548998132942739e-28
chr21   14181709        0.0013937899454307622   6.982887626608119e-07
chr21   14205523        0.00623506933299282     0.00014918022985855182
chr21   14208506        0.046374596221177       0.00028751565038632283
chr21   14214359        0.22458526252309008     0.001602013191933969
chr21   14215892        2.6458901468483425      0.005658162787052478
chr21   14216454        32.563740869346056      0.02395898515562496
chr21   14216846        2.296593902965297       0.024859249965587357
chr21   14216958        0.007906830684325574    0.024860135530624
chr21   14219368        0.0010110490694142274   0.02486257215888129

I don't know what I'm doing wrong, but the re_nodes.size always ends up being zero, which I understand as not having any recombination nodes. Any help would be appreciated!!

[petrelharp]

Hi, @jikhashkya - I think what's going on is that a massive chunk of the chr21 recombination map has zero recombination rate (in the chunk you gave, the first 14180791bp has no recombination). So, not only are there no RE events, but there are very few recombinations at all (I did get two distinct trees on the region, though, so there was some recombination).

Here's some code that does simulate from a heterogeneous rate map and does get RE events, so if there's a problem it's not general to rate maps.

rate_map = msprime.RateMap(position=[0, 50, 100], rate=[1, 2])
ts = msprime.sim_ancestry(samples=4, population_size=1, recombination_rate=rate_map, record_full_arg=True, random_seed=x)
ts.num_trees
# 100
np.sum(ts.tables.nodes.flags == msprime.NODE_IS_RE_EVENT) 
# 9590

Perhaps you could try a different chromosome with more recombination and see if that behaves as you expect?

[jikhashkya]

Will try that out. Do you think we would have observed more recombinations if the tool allowed a notion of "generations to simulate" as in other simulation tools ? It doesn't apply to coalescent simulators I suppose but I'm just curious.

[petrelharp]

Well, increasing the (effective) population_size will do just that, and I tried increasing it a lot, but still didn't get any. But yes, in general, this increases the number of recombination events.

[jeromekelleher]

Just to add to @petrelharp's answers here @jikhashkya - is there a particular reason for using a population_size of 1 here? If you're looking to model humans this is far too small and will given you several orders of magnitude less recombination than you would expect. A population_size of 10_000 is the usual rough number used, but much better estimates of human demography are available, e.g. in stdpopsim

[jikhashkya]

@jeromekelleher I think that was a mistake on my part where I understood it to be the number of different populations we would want to simulate. While stdpopsim does provide a better way to simulate human population, I want to also record the ARG for the simulation but I didn't see any clear options on doing that when using stdpopsim ? Would you happen to know if it's possible ?

[jeromekelleher]

Unfortunately not right now @jikhashkya, the released version of stdpopsim is based on msprime pre 1.0 and it's quite messy getting the information out.

I'm not sure running a full ARG simulation on a whole human chromosome is going to be feasible though in any case, it might result in a very large output.

[jeromekelleher]

The demography section of the documentation may be helpful here.