How to deal with recombination map length longer than chromosome length?

[lakishadavid]

My work requires chromosome-aware simulations, but I still would like to use the genetic maps managed through stdpopsim. My work focuses on generations 4 - 12 TMRCA. I'm attempting to create an extra-long chromosome that equals the total length of chromosomes 1 - 5 with four math.log(2) breaks representing where separations would be for five individual chromosomes. Because I need to convert the resulting tree sequence to VCF for further analysis, I'm creating a conversion dictionary with the appropriate positions so that the VCF will have the corresponding chromosome and position. The conversion dictionary uses rate_map.sequence_length. However, when I run rate_map = genetic_map.get_chromosome_map(f"chr{chromosome}"), I get the error message UserWarning: Recombination map has length 249218992.0, which is longer than chromosome length 248956422. The latter will be used. and I don't know how to account for that. Should I use the lengths based on the final positions in positions = [x + last_position + 1 for x in rate_map.position] instead of rate_map.sequence_length? How would I programmatically capture this warning message so that I can use the appropriate positions when converting the positions from the extra-long chromosome to the individual chromosomes?

species = stdpopsim.get_species("HomSap")
genetic_map = species.get_genetic_map("HapMapII_GRCh37")

starting_chromosome = 1
ending_chromosome = 5

r_break = math.log(2)
map_positions = []
map_rates = []
chrom_positions = [0]
genetic_map_positions = {} # as a check for the chromosome positions when converting to multiple chromosomes
conversion_dict = {} # used to convert the chromosome positions to the original positions in the VCF file

for chromosome in range(starting_chromosome, ending_chromosome + 1):
    
    # get_chromosome_map
    # Returns the genetic map for the chromosome with the specified ``id``
    # rtype: :class:`msprime.RateMap`
    # return: A :class:`msprime.RateMap` object
    rate_map = genetic_map.get_chromosome_map(f"chr{chromosome}")
    # genetic_map_positions is a check for the chromosome positions when converting to multiple chromosomes in vcf
    genetic_map_positions[chromosome] = int(rate_map.sequence_length)
    
    if chromosome != starting_chromosome:
        
        # map_positions
        # To make the simulation chromosome-aware, need to create one "long chromosome" of the total needed length
        # Get the last position in map_positions, which will be the last position in the previous map
        # Add the value of the last position + 1 to each value in the current map to make it sequential
        # Because the first position was 0, need the + 1 to avoid duplicate positions where maps join
        # Finally, add the positions to map_positions
        positions = [x + last_position + 1 for x in rate_map.position] # + 1 because it starts at zero, avoid duplicate numbers
        map_positions.extend(positions)
        conversion_dict[(int(last_position + 1), int(last_position + 1 + int(rate_map.sequence_length)))] = chromosome
        # last_position helps to build the sequential "long chromosome" of the total needed length of multiple chromosomes
        last_position = map_positions[-1]
        chrom_positions.append(int(last_position))
        
        
        # map_rates
        # Because each map starts with a span rate of nan, need to add the mean rate to the first span
        # To do so, convert the rate to a list, add the mean rate to position 0
        # And then convert it back to the array of the original shape
        # Finally, add the rates to map_rates
        shape = rate_map.rate.shape
        rates = rate_map.rate.tolist()
        rates[0] = rate_map.mean_rate
        rates = np.array(rates)
        rates.reshape(shape)
        map_rates.extend(rates)
        
    else:
        map_positions.extend(rate_map.position)
        last_position = map_positions[-1]
        chrom_positions.append(int(last_position))
        map_rates.extend(rate_map.rate)
        conversion_dict[(int(0), int(last_position))] = chromosome
    
    if chromosome != ending_chromosome:
        # set the recombination rate in the base pair segments separating chromosomes to log(2)
        # see https://tskit.dev/msprime/docs/latest/ancestry.html#example-2-simulating-without-a-pedigree
        map_rates.append(r_break)

rate_map = msprime.RateMap(position=map_positions, rate=map_rates)
assembly_len = int(rate_map.sequence_length)

print(f"The assembly length is {assembly_len:,}.")
print(chrom_positions)
print(genetic_map_positions)
print(conversion_dict)

Output (note that my conversion dictionary output uses the recombination map length instead of the chromosome length and so has the wrong positions):

/home/lakishadavid/.local/lib/python3.10/site-packages/stdpopsim/genetic_maps.py:127: UserWarning: Recombination map has length 249218992.0, which is longer than chromosome length 248956422. The latter will be used.
  warnings.warn(
/home/lakishadavid/.local/lib/python3.10/site-packages/stdpopsim/genetic_maps.py:127: UserWarning: Recombination map has length 243090997.0, which is longer than chromosome length 242193529. The latter will be used.
  warnings.warn(
/home/lakishadavid/.local/lib/python3.10/site-packages/stdpopsim/genetic_maps.py:127: UserWarning: Recombination map has length 191028645.0, which is longer than chromosome length 190214555. The latter will be used.
  warnings.warn(
The assembly length is 1,063,172,456.
[0, 249218992, 492309990, 690605550, 881634196, 1063172456]
{1: 249218992, 2: 243090997, 3: 198295559, 4: 191028645, 5: 181538259}
{(0, 249218992): 1, (249218993, 492309990): 2, (492309991, 690605550): 3, (690605551, 881634196): 4, (881634197, 1063172456): 5}

[andrewkern]

Hi @lakishadavid -- that warning about recombination map is being thrown because you are using

genetic_map = species.get_genetic_map("HapMapII_GRCh37")

but the default human genome is stdpopsim is GRCh38 currently. If you swap out that map for the h38 map that warning will go away, e.g.,

genetic_map = species.get_genetic_map("HapMapII_GRCh38")

Let me know if that fixes it for you?

cheers,
Andy

[lakishadavid]

That fixed it. Thanks!

[andrewkern]

Great!