Reaching equilibrium with recapitation

[Aude Caizergues]

Hello,

I have a question regarding the recapitation step in my simulations, I know it is not really a SLiM question but I don't see where else to ask the question...

I want to run a series of non-WF model with different parameters of life history traits and environment to investigate how each would affect patterns of neutral evolution. To do so, I run my simulation in SLiM and then overlay the neutral mutations with the recapitation process with pyslim.

At the beginning of my project, I was running my simulations with 5,000,000 generations burn-in, which is really time consuming, so I decided to re-run my analyses with only 10,000 burn-in. So I ended up with two similar models differing just in the number of generation burn-in before I started the tree recording (i.e. I am recording 5100 generations before the end of my simulations). The models are strictly similar but when in move on the the recapitation process, I use the same parameters but still end up with different levels of diversity in my two scenarios: the 5 million generation burn-in has about 115,000 mutations whereas the 10,000gen burn-in has about 95,000 (I ran some replicates). And for the 10,000gen, the diversity keeps on decreasing with time whereas it is constant in the 5million gen, see the graphs below (please don't mind the colors, both lines are just different subsets of individuals):

I don't understand where this behaviour comes from as the previous burn-in history is not recorded in the tree. I'm I missing something ?

Thank you in advance for your help,

Aude

[Ben Haller]

Hi Aude.  This question is probably one for Peter, not me, but perhaps I can guess what his first question will be: could you please post the scripts for these two alternative approaches, so that we can see exactly what you are doing?  :->

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University


Aude Caizergues wrote on 2/6/24 10:39 AM:

--
SLiM forward genetic simulation: http://messerlab.org/slim/
---
You received this message because you are subscribed to the Google Groups "slim-discuss" group.
To unsubscribe from this group and stop receiving emails from it, send an email to slim-discuss...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/slim-discuss/a1a64d99-1442-4ab5-8d9c-d0b3ad5cb1c0n%40googlegroups.com.

[Aude Caizergues]

Yes sure ! The slim scripts are attached and this is the script I used for the recapitation:

genmean = np.nanmean(pyslim.individual_ages_at(orig_ts,0))
recap_ts = pyslim.recapitate(orig_ts,recombination_rate=1e-9, ancestral_Ne=5000, random_seed=42)
ts = msprime.sim_mutations(
    recap_ts,
    rate=(1e-8)/genmean,
    model=msprime.SLiMMutationModel(type=0),
    keep=True)
 
sum([t.num_roots == 1 for t in ts.trees()])
print(f"The tree sequence now has {ts.num_mutations} mutations,\n"
      f"and mean pairwise nucleotide diversity is {ts.diversity():0.3e}.")

I must also add, that each time I have been trying to recapitate a model with 10,000 I've observed decreasing pi over time. So I though it could come from a wrong Ne input. To have an idea of the actual Ne,  I ran simulations with neutral mutations in SLiM and fount a Ne of about 800. But even with recapitating with Ne=800 I obtained a decreasing pi with time (even if the pi was already lower than what I observed with the 500,000 burin-in (and the 500,000 had a stable-ish pi)). So I guess there is a bug somewhere that I can't see...

Thank you !

PS: I said 5,000,000 burn-in in my previous message, it is actually 500,000.

[Peter Ralph]

Hello, Aude -

Well, I think you're thinking in the same direction I am here: if the recapitation portion creates initial diversity at the start of the burn-in period that is lower than it should be, then you'd see a pattern like what you're seeing. I note that Ne is about 10K, and in the plots you show, the decrease happens at about 5K generations ago, which is (I think) 10K generations after the start of burn-in, which is reasonable for when the effect of too-low diversity from recapitation would show up. This is easy to check: just run the demographic model that you are using for recapitation all by itself and compute diversity, and compare that to diversity in the very long simulation.

As for why you might not be having enough diversity produced by recapitation: there are tricky things to do with time scales. I see you've incorporated generation time,

but I don't think it's quite right? At least, it's not done in the same way as we outline here:

  https://tskit.dev/pyslim/docs/latest/time_units.html#example-calculating-generation-times-in-a-slim-simulation

(in particular, you haven't scaled Ne).

best of luck,

  Peter

p.s. You said that "I am recording 5100 generations before the end of my simulations" and  "the previous burn-in history is not recorded in the tree." But, unless I'm missing something, tree sequence recording is actually on for the entire simulation (i.e., including burn-in), so the trees are recorded; just you are not Remembering individuals for that time period.

---

From: slim-d...@googlegroups.com <slim-d...@googlegroups.com> on behalf of Aude Caizergues <audeemilie...@gmail.com>
Sent: Tuesday, February 6, 2024 9:01 AM
To: slim-discuss <slim-d...@googlegroups.com>
Subject: Re: Reaching equilibrium with recapitation

 

To view this discussion on the web visit https://groups.google.com/d/msgid/slim-discuss/a86b5c36-440f-460a-842d-9bf3ed567730n%40googlegroups.com.

[Aude Caizergues]

Hello,

My apologies for the late update. 

You were right about your guess, I indeed completely missed out on the scaling of Ne... I computed the generation time following the tutorial in the link provided and got a rough estimate of Ne with a long simulation including mutations and then computing pi and scaled my Ne from that. It worked ! I was able to generate the good amount of diversity and reach equilibrium. 

Thank you for your help,

Best regards,

Aude