Stdpopsim: Tennessen model QC

[ndukler]

Just opening this up as a forum for the QC of the Tennessen two population OOA model as per the earlier discussion.

[jeromekelleher]

I've just assigned @ndukler and @carjed as the people directly involved in this. Is this correct?

[ndukler]

Yep, thanks! I'll also migrate my previous comments about QC to this thread.

[ndukler]

So I've been writing the QC script for this model, and it kinda quickly became obvious that, for whatever reason, working from the primary paper, the model is somewhat ambiguous. There are multiple parameter sets, multiple possible interpretations of which migration parameters should be included etc. For example, it's not obvious which of the three parameter sets inferred in Gavel was used in Tennessen (although it seems to be the low-coverage + exons), nor is it obvious what happens to the migration rates with respect to the first and second population growth periods. I'll bring this up at the next chat, just thought I'd make a note of it here in case anyone had any ideas.

[ndukler]

Additionally, my calculations lead to a somewhat different current European population size than what Tennessen shows in figure 2B. It is entirely possible that I am just bad at this, but the African calculation seems to be correct so it makes me a little suspicious.

[ndukler]

As per @carjed the population ID numbers are AF=0 and EU=1, to be consistent with the GutenkunstThreePop model, which uses AF=0, EU=1, and AS=2.

[ndukler]

@jeromekelleher As I kinda expected, my implementation of the model did not match the production version. Would it be a good idea to have a diff() function that would report model differences per epoch in a format similar to the DemographyDebugger() program? Otherwise I'm just comparing two DemographyDebugger() outputs which is doable, but feels clumsy.

[carjed]

I agree that we should assume the Tennessen paper uses the low-coverage + exons parameter set from Gravel et al., because that is the only one that lines up with the OOA migration at 51kya (and this is the best-fitting model reported in Fig. 4 of Gravel et al.).

The Tennessen supplement states that they fixed the more ancient parameters of the Gravel model, "with three added parameters to allow for a recent growth, namely a time for growth change and two recent growth rates (one in the European and one in the African population)", so all other parameters in the low-coverage + exons column in Table 2 of Gravel should be fixed, giving us migration parameters of 15e-5 and 2.5e-5 before and after the first European bottleneck, respectively.

Looking closer at the debugging output of my implementation (see below), I am also getting slightly different values of Ne for both AFR and EUR (~432,000 and ~501,000, respectively, compared to 424,000 and 512,000 in Fig. 2B of Tennessen), but the math works out to get the appropriate Ne values at 5.1kya and 23kya.

For example, prior to the recent growth starting at 5.1kya, Ne_AFR = 14474, so the present Ne for AFR should be:

>>> 14474/math.exp(-0.0166*5115/25)
432124

...so I'm really not sure where they are getting 424,000, unless they A) used a lower Ne_AFR (in the ballpark of 14,200), B) used a growth rate slower than 1.66%, or C) the time at which growth started was more recent than 5115ya. I've played around with different combinations of rounding errors, potential typos, etc., and the closes I can come to recapitulating an Ne of 424,000 is setting a growth rate of 1.65%, which gives an Ne of 423,373.

>>> TennessenTwoPop_model.debug()
=============================
Epoch: 0 -- 204.6 generations
=============================
     start     end      growth_rate |     0        1
   -------- --------       -------- | -------- --------
0 |4.32e+05 1.45e+04         0.0166 |     0        0
1 |5.01e+05 9.28e+03         0.0195 |     0        0

Events @ generation 204.6
   - Migration rate change for (0, 1) to 2.5e-05
   - Migration rate change for (1, 0) to 2.5e-05
   - Population parameter change for 1: initial_size -> 9279.212349452768 growth_rate -> 0.00307
   - Population parameter change for 0: initial_size -> 14474 growth_rate -> 0
=================================
Epoch: 204.6 -- 920.0 generations
=================================
     start     end      growth_rate |     0        1
   -------- --------       -------- | -------- --------
0 |1.45e+04 1.45e+04              0 |     0     2.5e-05
1 |9.28e+03 1.03e+03        0.00307 |  2.5e-05     0

Events @ generation 920.0
   - Migration rate change for (0, 1) to 0.00015
   - Migration rate change for (1, 0) to 0.00015
   - Population parameter change for 1: initial_size -> 1861 growth_rate -> 0
==================================
Epoch: 920.0 -- 2040.0 generations
==================================
     start     end      growth_rate |     0        1
   -------- --------       -------- | -------- --------
0 |1.45e+04 1.45e+04              0 |     0     0.00015
1 |1.86e+03 1.86e+03              0 |  0.00015     0

Events @ generation 2040.0
   - Mass migration: lineages move from 1 to 0 with probability 1.0
===================================
Epoch: 2040.0 -- 5920.0 generations
===================================
     start     end      growth_rate |     0        1
   -------- --------       -------- | -------- --------
0 |1.45e+04 1.45e+04              0 |     0     0.00015
1 |1.86e+03 1.86e+03              0 |  0.00015     0

Events @ generation 5920.0
   - Population parameter change for 0: initial_size -> 7310
================================
Epoch: 5920.0 -- inf generations
================================
     start     end      growth_rate |     0        1
   -------- --------       -------- | -------- --------
0 |7.31e+03 7.31e+03              0 |     0     0.00015
1 |1.86e+03 1.86e+03              0 |  0.00015     0

[carjed]

Also noticed I made a stupid error in my previous PR (#41)—the latest commit in my fork, https://github.com/carjed/stdpopsim/commit/a9eea76f150dd125195ce14649be7929db509780, resolves this, and the debugging output above is based on this revision.

[ndukler]

For my modern population size estimates I got:
-N_AF ~ 432,124
-N_EU ~ 501,436
which agree with your estimates but disagree with the original paper.

[ndukler]

Here's my debugger output. The main point on which we differ is the modern migration setting I think:

=============================
Epoch: 0 -- 204.6 generations
=============================
     start     end      growth_rate |     0        1    
   -------- --------       -------- | -------- -------- 
0 |4.32e+05 1.45e+04         0.0166 |     0     2.5e-05 
1 |5.01e+05 9.28e+03         0.0195 |  2.5e-05     0    

Events @ generation 204.6
   - Population parameter change for 0: initial_size -> 14474 growth_rate -> 0 
   - Population parameter change for 1: initial_size -> 9279.212349452768 growth_rate -> 0.00307 
   - Migration rate change for (0, 1) to 2.5e-05
   - Migration rate change for (1, 0) to 2.5e-05
=================================
Epoch: 204.6 -- 920.0 generations
=================================
     start     end      growth_rate |     0        1    
   -------- --------       -------- | -------- -------- 
0 |1.45e+04 1.45e+04              0 |     0     2.5e-05 
1 |9.28e+03 1.03e+03        0.00307 |  2.5e-05     0    

Events @ generation 920.0
   - Population parameter change for 1: initial_size -> 1032 growth_rate -> 0 
   - Migration rate change for (0, 1) to 0.00015
   - Migration rate change for (1, 0) to 0.00015
==================================
Epoch: 920.0 -- 2040.0 generations
==================================
     start     end      growth_rate |     0        1    
   -------- --------       -------- | -------- -------- 
0 |1.45e+04 1.45e+04              0 |     0     0.00015 
1 |1.03e+03 1.03e+03              0 |  0.00015     0    

Events @ generation 2040.0
   - Mass migration: lineages move from 1 to 0 with probability 1.0
===================================
Epoch: 2040.0 -- 5920.0 generations
===================================
     start     end      growth_rate |     0        1    
   -------- --------       -------- | -------- -------- 
0 |1.45e+04 1.45e+04              0 |     0     0.00015 
1 |1.03e+03 1.03e+03              0 |  0.00015     0    

Events @ generation 5920.0
   - Population parameter change for 0: initial_size -> 7310 
================================
Epoch: 5920.0 -- inf generations
================================
     start     end      growth_rate |     0        1    
   -------- --------       -------- | -------- -------- 
0 |7.31e+03 7.31e+03              0 |     0     0.00015 
1 |1.03e+03 1.03e+03              0 |  0.00015     0

[carjed]

Yeah, Fig. S5 in the Fu paper is somewhat ambiguous about the AFR<-->EUR migration rates from 0-5.1kya so I set it to 0, but the Gravel model implies that the migration rate of 2.5e-5 carries through to the present, so fixing this makes sense to me.

Another difference appears to the European bottleneck from Ne=1,861 -> 1,032 @ generation 920.0. Yours seems to be carrying through Ne=1,032 until the Europeans merge back into the ancestral population.

[jeromekelleher]

@jeromekelleher As I kinda expected, my implementation of the model did not match the production version. Would it be a good idea to have a diff() function that would report model differences per epoch in a format similar to the DemographyDebugger() program? Otherwise I'm just comparing two DemographyDebugger() outputs which is doable, but feels clumsy.

I've copied this to a new issue #54 and started discussion there.

[ndukler]

Yeah, Fig. S5 in the Fu paper is somewhat ambiguous about the AFR<-->EUR migration rates from 0-5.1kya so I set it to 0, but the Gravel model implies that the migration rate of 2.5e-5 carries through to the present, so fixing this makes sense to me.

Another difference appears to the European bottleneck from Ne=1,861 -> 1,032 @ generation 920.0. Yours seems to be carrying through Ne=1,032 until the Europeans merge back into the ancestral population.

Yup, you're right, the European bottleneck numbers were a bug on my part. The lack of an additional population for additional people to appear from threw me off, thanks!

[ndukler]

So where do we go from here? I can contact Joshua Akey to check the migration thing, but even after that the current setup makes it hard to run verify() since the candidate production model is always on someone else's fork unless it has been pulled into master which seems premature without QC. The two possible solutions are the I just clone @carjed 's fork, or we have a "pre-production" branch on the central repository onto which candidate models can be pulled -> QC'ed -> Merged into master.

[jeromekelleher]

A simpler expedient might be to just not document/offer CLI support for any models that haven't been QC'd fully yet. Maintaining multiple branches is a major PITA and it always gets messy. I've opened an issue in #55 to start formalising the QC process.

[ndukler]

@jeromekelleher How should I add the QC model to the unit tests? I ask b/c I wrote my QC more as a notebook with detailed provenance for the parameters. Should I just wrap the whole thing in a function called test_demographic_model and add it to the appropriate class in tests/test_homo_sapiens.py ? The only downside to that is that is will add a lot of code to the test file.

[carjed]

@ndukler, here's what I proposed to keep things tidy—basically we'll create a homo_sapiens_qc.py module in the qc directory containing the blind implementation, import this module into tests/test_homo_sapiens.py, and add the verify_equal unit test: https://github.com/popgensims/stdpopsim/issues/55#issuecomment-481889395

This way we're keeping the production models (stdpopsim/{species}.py), the QC models (qc/{species}_qc.py), and the unit tests (tests/test_{species}.py) in their own separate silos.

[jeromekelleher]

I think @carjed's plan is an excellent one.

[ndukler]

Ok, so I've implemented QC as a module and ported my implementation over to that format as discussed in #55. I haven't added the unit test yet b/c it would fail right now due differences in the migration matrix. I emailed Josh and he referred me to a post-doc of his who works heavily with the Tennessen model so I'm waiting for them. When I get that back I'll squash the commits and start a pull request.

[carjed]

@ndukler, maybe go ahead and open a PR, even though the unit test will fail? This way we can try out the proposed workflow (https://github.com/popgensims/stdpopsim/issues/55#issuecomment-482217521) for implementing hotfixes in the production model.