Hello,
Before testing an empirical model in PSMC, we are trying to simulate a simple demographic model to make sure we are getting out of PSMC what we believe should match the simulated model. One of the main problems we're having now is getting the x-axis time value inferred from the psmc to match what we simulate in msprime. For example, we have this model:
demographic_model = msprime.Demography()
demographic_events=[
demographic_model.add_population(initial_size = 5000, growth_rate = 0), #Size at t=0, then shrinks to ~11k individuals by gen. 40k
demographic_model.add_population_parameters_change(time=600000, growth_rate = 0, population = 0), #Stop growth at 200k gens ago
demographic_model.add_population_parameters_change(time=1000000, initial_size=10000, population = 0)
]
Output from the demography debugger:
I assume this means the population's origin should go back to 1e6 or 6e5 generations ago. We use these values:
mut_rate = 8e-8
r_chrom = 1e-7
chr_length = int(1e7)
But the psmc output looks like this:
It's a little curious that the Ne is ~10x less than that simulated, and the x-axis is ~10x less than what it should be. But I'm not quite sure why. This model used multiple (independently recombining) chromosomes, but we get similar results with a model simulating a single chromosome. I wonder if there aren't enough recombination events to infer the time that long ago? I set the chromosome length and rec. rate to be the inverse so there is ~1 recom/chr. expected.
What am I missing? Any insight would be greatly appreciated!
Thanks in advance,
Jared