Multiple chromosomes with pedigree builder

[lakishadavid]

I'm successfully using msprime's pedigree-aware simulations if I use chromosome 1. However, simulating first cousins (specifically, first cousins to eighth cousins in a 10-generation pedigree) is important for my project. The expected amount of sharing for second cousins is on average 212.50 cM, so I could make do with just simulating with chromosome 1, but this does not reflect my real-life sampling plan of sampling first cousins to reduce uncertainty in a 10-generation inferred pedigree. The expected amount of sharing for first cousins is on average 850.00 cM so I would need to simulate multiple chromosomes. I run into errors as I try to do so and would like your advice on how to proceed.

UPDATE: I now have this part of the question working. See #2081 (comment). I'm now only trying to figure out the second part of this question about the missing regions error.

First, I confirmed that my code works using chromosome 1 and the NIH HapMap 37 genetic map. All is well from building the pedigree to running the simulation, ibd_segments, and writing a VCF file. (Pedigree building includes assigning specific populations among AFR, ADMIX, or EUR to specific individuals as per the pedigree definition.) Then, I followed your example here adapting it to simulate five chromosomes. This runs but when I used keep_intervals post-simulation based on the intervals in chrom_positions, only the ones for the first interval were present. I get errors on the subsequent intervals stating that they did not exist. The VCF file produced from this attempt only had positions for the first interval as well.

import math

r_chrom = 1e-8
r_break = math.log(2)
# from https://www.ncbi.nlm.nih.gov/grc/human/data?asm=GRCh37
chrom_positions = [0, 249250621, 492449994, 690472424, 881626700, 1062541960]
map_positions = [
    chrom_positions[0],
    chrom_positions[1],
    chrom_positions[1] + 1,
    chrom_positions[2],
    chrom_positions[2] + 1,
    chrom_positions[3],
    chrom_positions[3] + 1,
    chrom_positions[4],
    chrom_positions[4] + 1,
    chrom_positions[5]
]
rates = [r_chrom, r_break, r_chrom, r_break, r_chrom, r_break, r_chrom, r_break, r_chrom]
rate_map = msprime.RateMap(position=map_positions, rate=rates)

Missing regions error: I also tried to use this strategy with the genetic maps by renumbering the positions gathered from the map as if I'm working with one extra-long chromosome pre-simulation, and then selecting the appropriate intervals using keep_intervals and renumbering them back post-simulation according to the specific chromosome that those positions represent pre-simulation. I don't know if it makes sense to do this renumbering scheme. In any event, when I try to run the simulation, I'm getting the error # Missing regions of the genome other than the flanks are currently not supported..

import math

r_break = math.log(2)
map_positions = []
map_rates = []

for chromosome in range(1,6):
    map_file_name = os.path.join(directory, f"starbright/references/nih_hapmap37_map/genetic_map_GRCh37_chr{chromosome}.txt")
    rate_map = msprime.RateMap.read_hapmap(map_file_name)
    
    if chromosome != 1:
        last_position = map_positions[-1]
        positions = [x + last_position + 1 for x in rate_map.position] # + 1 because it starts at zero, avoid duplicate numbers
        map_positions.extend(positions)
    else:
        map_positions.extend(rate_map.position)

    map_rates.extend(rate_map.rate)
    if chromosome != 5:
        map_rates.append(r_break)

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

Could you help me to understand what I'm not getting quite right?

[apragsdale]

Hi @lakishadavid -- I'm not really sure what's going wrong with setting up the multi-chromosome map in the first code block. Nothing pops out to me as wrong. However, since you were able to get the simulation working with the FixedPedigree model, you might be able to altogether avoid this whole process of stitching together the recombination maps.

With independent simulations through the same fixed pedigree, those simulations will ensure that you have free recombination between those simulated chromosomes. You would only need to stitch together multiple chromosomes into a single map if your simulations are not conditioned on a pedigree (if you used the dtwf model instead, for example). With the luxury of a fixed pedigree, that whole recombination-map-stitching process is unnecessary.

So I would suggest running each chromosome with its own recombination map independently under the same input pedigree. You might need to be careful about how the simulation is finished off above the pedigree, if you need trees to be fully coalesced. But parallelizing the simulations like this should come with the added bonus of being more efficient and simplifying your code. I hope this is helpful!

[apragsdale]

(Also, now that the fixed pedigree model is in the main release, we should update those multiple chromosome docs to explain this as a cleaner option for simulating multiple chromosomes.)

[lakishadavid]

I don't see how running independent simulations for each chromosome would not produce more IBD sharing than expected when I then go to concatenate the VCF files. When or how are you suggesting that I run the independent simulations through the same fixed pedigree such that this issue is avoided?

Here is how I finish the pedigree

# Build the msprime pedigree

# Set so that the follow-up recapitation simulations will use the defined 
# demography model. Necessary when drawing founders from multiple populations
pb = msprime.PedigreeBuilder(
    demography,
    individuals_metadata_schema=tskit.MetadataSchema.permissive_json()
    )

# Build out the pedigree using input text pedigree
pb, txt_ped_to_tskit_key = add_individuals_to_pedigree(pb, text_pedigree)

# Set the initial state of tree sequence
# option 1, set sequence_length = len
pedigree = pb.finalise(sequence_length = assembly_len)

And here is how I run the simulation.

t1_start = process_time()

# simulation within fixed pedigree
ped_ts = msprime.sim_ancestry(
    initial_state = pedigree,
    recombination_rate = rate_map,
    model = "fixed_pedigree",
    record_provenance=False,
    random_seed = seed + 100
    )

# complete the simulation to simulate ancestry older than what is defined in the pedigree
# and because drawing from multiple populations 
completed_ts = msprime.sim_ancestry(
    initial_state = ped_ts,
    recombination_rate = rate_map,
    demography = demography,
    model=[
        msprime.DiscreteTimeWrightFisher(duration=20),
        msprime.StandardCoalescent(),
    ],
    record_provenance=False,
    random_seed = seed + 200
    )

# add mutations
ts = msprime.sim_mutations(
          completed_ts,
          rate = mutation_rate,
          random_seed = seed + 300
          )

t1_stop = process_time()

runtime = t1_stop - t1_start
hours = runtime // 3600
minutes = (runtime % 3600) // 60
seconds = (runtime % 3600) % 60

print(f"Runtime: {hours:.0f}:{minutes:.0f}:{seconds:.0f}\n")

ts

[apragsdale]

I don't see how running independent simulations for each chromosome would not produce more IBD sharing than expected when I then go to concatenate the VCF files. When or how are you suggesting that I run the independent simulations through the same fixed pedigree such that this issue is avoided?

Sorry, I don't follow the question. Are you concerned that there will be more IBD than expected, or less? How does concatenating the VCF change things here?

If your single chromosome simulations (for just chromosome 1) are giving the correct amount of IBD, I don't see why simulating multiple chromosomes through the same pedigree would make the simulated IBD deviate from expectations.

For finishing off your simulations, if your pedigree extends back into the past long enough (maybe >10 generations?) switching to DTWF/Hudson should have a very minor effect on errors in correlations between chromosomes, as those should be well-captured by the pedigree.

[lakishadavid]

Let's say for a pair of second cousins: if I run this on one chromosome and get IBD sharing of 244 cM and then independently run it again on another chromosome and get IBD sharing of 180 cM, wouldn't the concatenated VCF file show 424 cM of IBD sharing? And even if so, there's a 20% probability that second cousins share 424 cM using real data so that's okay if I was only using two chromosomes, I guess, but I'm using five chromosomes so that I can also simulate first cousins. It seems that the cumulative IBD sharing would then be too much from independent runs for second cousins and that the simulation needs to be chromosome-aware. I know msprime does not have the chromosome-aware functionality, but conceptually, independent runs do have a meaningful bearing on the totals.

Even if we disagree on this concern, none of this explains why the simplified example for multiple chromosomes from the manual does not produce segments beyond the first interval when investigating it with keep_intervals.

[apragsdale]

I think I'm not understanding this. Expected total IBD sharing will be a function of relatedness as well as genome length, and the number and size of each block will depend on the chromosome/recombination structure. But in any case, simulating each chromosome independently through the pedigree ensures free recombination, which is what you want for multi-chromosome simulations.

I tested out the code you pasted to generate the initial recombination map, and then parsed out each chromosome as in the example in the docs. I didn't run into any errors and got the list of five tree sequences with the correct lengths in the end, so I'm not really sure what's going on there.

[lakishadavid]

Update: I now have the multiple chromosomes from the adapted example working. I was missing the trim() method from keep_intervals. This produces the appropriate sequence length for chromosome 5, etc.

ts5 = ts.keep_intervals([[chrom_positions[4], chrom_positions[5]]], simplify=False).trim()
ts5

I'm now left to figure out how to correct the # Missing regions of the genome other than the flanks are currently not supported. error from the code in my original question when I try to use the HapMap genetic maps instead.

[jeromekelleher]

Which code chunk is throwing this error @lakishadavid? I think what's happening is that we're merging several maps with regions of unknown recombination rate at the flanks, which is giving a single map lit lots of unknown rates dotted through it. This is discussed in #1604.

Probably the best way to get this working would be to update your maps to insert 0 recombination rates explicitly at the start and end. However you need to either manually do delete intervals at the end, or be very careful of interpreting results from these large zero recombination rate chunks.

[lakishadavid]

Using np.argwhere(np.isnan(rate_map.rate)), the first position of each map was nan (as you mentioned). Following, @hyanwong's suggestion at #1604 (comment), I replaced all nan for chromosomes 2 - 5 with rate_map.mean_rate. (Thank you for having these open discussions.) It now runs fine, and I can successfully subset the resultant tree sequence on interval ranges.

Here is the full code for stitching together chromosomes 1 - 5 (with chrom_positions used later in keep_intervals in preparation for converting the positions back to the original positions in the VCF file):

import math

r_break = math.log(2)
map_positions = []
map_rates = []
chrom_positions = [0]

for chromosome in range(1,6):
    map_file_name = os.path.join(directory, f"starbright/references/nih_hapmap37_map/genetic_map_GRCh37_chr{chromosome}.txt")
    rate_map = msprime.RateMap.read_hapmap(map_file_name)
    
    if chromosome != 1:
        last_position = map_positions[-1]
        chrom_positions.append(last_position)
        positions = [x + last_position + 1 for x in rate_map.position] # + 1 because it starts at zero, avoid duplicate numbers
        map_positions.extend(positions)

        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)
        map_rates.extend(rate_map.rate)
    
    if chromosome != 5:
        map_rates.append(r_break)

rate_map = msprime.RateMap(position=map_positions, rate=map_rates)
chrom_positions.append(rate_map.sequence_length)

[jeromekelleher]

Great!