search

popsim-consortium / stdpopsim

A library of standard population genetic models
GNU General Public License v3.0
125 stars 87 forks source link

`msprime.RateMap.slice` throws error when contig covers portion of genetic map with no data #1387

Open nspope opened 2 years ago

nspope commented 2 years ago

For example:

$ stdpopsim HomSap -g HapMapII_GRCh38 -c chr22 -o foo.ts --right 10000000 -d OutOfAfrica_2T12 AFR:2 EUR:3
Traceback (most recent call last):
  File "/home/natep/.miniconda3/envs/stdpopsim/bin/stdpopsim", line 8, in <module>
    sys.exit(stdpopsim_main())
  File "/home/natep/projects/stdpopsim/stdpopsim/stdpopsim/cli.py", line 1136, in stdpopsim_main
    run(args)
  File "/home/natep/projects/stdpopsim/stdpopsim/stdpopsim/cli.py", line 1127, in run 
    args.runner(args)
  File "/home/natep/projects/stdpopsim/stdpopsim/stdpopsim/cli.py", line 749, in run_simulation
    contig = species.get_contig(
  File "/home/natep/projects/stdpopsim/stdpopsim/stdpopsim/species.py", line 229, in get_contig
    return stdpopsim.Contig.species_contig(
  File "/home/natep/projects/stdpopsim/stdpopsim/stdpopsim/genomes.py", line 506, in species_contig
    recomb_map = recomb_map.slice(left=left, right=right, trim=True)
  File "/home/natep/.miniconda3/envs/stdpopsim/lib/python3.10/site-packages/msprime/intervals.py", line 428, in slice
    return RateMap(position=position - left, rate=rate)
  File "/home/natep/.miniconda3/envs/stdpopsim/lib/python3.10/site-packages/msprime/intervals.py", line 86, in __init__
    raise ValueError("All intervals are missing data")
ValueError: All intervals are missing data

which happens because there's no recombination rates for the first 15 Mb in the chr22 map.

nspope commented 2 years ago

RateMap.slice will also fail if the genetic map is shorter than the chromosome (if the chromosome is shorter than the genetic map, the latter will be clipped). I don't think that's wrong, necessarily, but the error message may be a bit obtuse.