Telomere region is not removed in SLiM simulation

[daikitag]

Genetic map in stdpopsim is read by using read_hapmap function in msprime (https://github.com/popsim-consortium/stdpopsim/blob/main/stdpopsim/genetic_maps.py#L105), and is loaded to SLiM by using msprime_rm_to_slim_rm function, (https://github.com/popsim-consortium/stdpopsim/blob/main/stdpopsim/slim_engine.py#L884). The genetic map loaded by msprime has NaN at the beginning to account for telomere regions, but it is converted into a recombination rate of 0 in SLiM, (https://github.com/popsim-consortium/stdpopsim/blob/main/stdpopsim/slim_engine.py#L908).

This long non-recombining region at the beginning gets removed in the final output from recapitation and simplifacation in msprime simulation, but it causes some issues when there is selection implemented in SLiM. Specifically, if users follow "Simulating with a genome-wide DFE" in stdpopsim manual (https://popsim-consortium.github.io/stdpopsim-docs/stable/tutorial.html#simulating-with-a-genome-wide-dfe), initializeGenomicElement will start from 0 (https://github.com/popsim-consortium/stdpopsim/blob/main/stdpopsim/slim_engine.py#L1327), and there is a long 0 recombination rate region at the beginning of the SLiM tree sequence output with a bunch of mutations. This will influence the final tree sequence output, even though there should not be any selected mutations in the telomere region.

Suggested fixes:
Instead of simply setting a recombination rate of 0 for NaN regions, we might be able to store the region in an array and remove that region from the DFE region.

[petrelharp]

Oh, good catch. Thanks for the nice report, and yes, this should be fixed. Hm: mind providing a short script showing this happening? (I know you referred to the code above, but a self-contained example showing there's mutations where there shouldn't be would be good.)

Let's see: as you say, genomic elements are initialized here, which simply loops over the DFEs.

Seems like the fix should come upstream: when we do contig.add_dfe( ), for instance, it could remove the missing intervals; and perhaps that would fix it? We'd also need to fix this here.

[nspope]

This came up before -- so previous conversation is maybe relevant? #1422 (comment)

[nspope]

The followup comment to the one I linked reminded me that we have a test checking that there's no selected mutations simulated in the missing flanks -- see #1617 -- are you seeing output that suggests this isn't working as anticipated, @daikitag ? (It could totally be the case that we'd thought this was fixed, when it isn't -- in which case the test should get fixed, too)

[daikitag]

@nspope @petrelharp
Thank you very much for your reply. I have attached the screenshot of the Jupyter Notebook, which clearly shows that there are mutations in the telomere region. I also checked the trees from the SLiM output in tsbrowse, and there were loads of mutations in telomere region.

I think the SLiM codes that were generated here, https://github.com/popsim-consortium/stdpopsim/blob/main/tests/test_slim_engine.py#L1809, don't have mutations outside of [left, right] region, because in that setting, the initializeGenomicElement starts from the left position instead of 0.

[daikitag]

This is the code that I used:

homsap = stdpopsim.get_species("HomSap")
contig = homsap.get_contig("chr22", genetic_map="HapMapII_GRCh38")
dfe = homsap.get_dfe("Gamma_K17")
contig.add_dfe(np.array([[0, contig.length]], dtype="int"), dfe)
dfe_breaks, dfe_type = contig.dfe_breakpoints()

engine = stdpopsim.get_engine("slim")
demogr = stdpopsim.PiecewiseConstantSize(1e3)

ts = engine.simulate(
    demographic_model=demogr,
    contig=contig,
    samples={"pop_0": 10},
    seed=1234,
    _recap_and_rescale=False,
)

[nspope]

ah, got it -- the previous fix (and test) are testing that when we subset a contig via left=... and right=... the DFE gets clipped, but this doesn't help when the recombination map for the entire chromosome has missing portions. Thanks for spotting this, and for the MWE!

So one simple way to fix this is by having Species.get_contig set left, right to the first/last nonmissing positions of the recombination map, I think?

[petrelharp]

Let's see: species.get_contig passes this to Contig.species_contig. That seems like it would work, if the only missing areas are always at the start and/or end of the contig? Alternatively, we could use the exclusion_mask? Seems like left and right are the way to go, though, if that works, because those aren't properties of the resulting Contig, as oppsed to exclusion_mask.

[nspope]

The exclusion mask is applied after simulation, I thought?

stdpopsim/stdpopsim/slim_engine.py

Lines 1749 to 1750 in e23d255

	if contig.exclusion_mask is not None:
	ts = stdpopsim.utils.mask_tree_sequence(ts, contig.exclusion_mask, True)

So wouldn't deal with the issue of unintended, selected mutations distorting nearby genealogies. And presumably we'd error out if there are NaN's in the interior of the ratemap -- this is what happens with msprime.

I'll go ahead and put together a PR to do the left/right fix, we can discuss more there