Stdpopsim: QC for PapuansOutOfAfrica_10J19 model

[andrewkern]

opening this issue from merge of https://github.com/popsim-consortium/stdpopsim/pull/372

Jacobs et al doi is https://doi.org/10.1016/j.cell.2019.02.035

[grahamgower]

Thanks @andrewkern.

The population indices for QCing can be found here:
https://stdpopsim.readthedocs.io/en/latest/catalog.html#sec_catalog_homsap_models_papuansoutofafrica_10j19

Most model parameters are taken from Jacobs et al., Table S5 and Figure S5, but a few were found in the supplement text (section S10.i and S10.h, pgs. 30/45 & 31/45).

[grahamgower]

I occurs to me that maybe no one wants to QC this because its big and will be time consuming. So I'll offer an incentive: if you QC this model, I promise to QC a future model of your choice.

[ndukler]

I'll do it this weekend. I've just been a bit bogged down with my other work recently.

[noscode]

I have almost done it. I have some questions though:
1) Time of Neanderthal Altai and Neanderthal Introgressing split event - I found in the Table S5A 3375 generations, but in implemented model there is 2058. It seems to be sampling time of DenA population.
2)

# Fig. 4B, and S10.h p. 30/45
T_Den1_Papuan_mig = 29.8e3 / generation_time
T_Den2_Papuan_mig = 45.7e3 / generation_time

I have not found that parameters but I think that they should be both equal to 1353. In section S10.h I found the following:

1. Instead of a single Denisovan introgression into Australians 1353 generations ago of 4%, there are two introgressions 1353 generations ago into Papuans, one from D1 and one from D2

[grahamgower]

I have almost done it. I have some questions though:

Awesome, thanks @noscode!

1. Time of Neanderthal Altai and Neanderthal Introgressing split event - I found in the Table S5A 3375 generations, but in implemented model there is 2058. It seems to be sampling time of DenA population.

Yep, that looks wrong. I'll send a pull request changing that.

# Fig. 4B, and S10.h p. 30/45
T_Den1_Papuan_mig = 29.8e3 / generation_time
T_Den2_Papuan_mig = 45.7e3 / generation_time

I have not found that parameters but I think that they should be both equal to 1353. In section S10.h I found the following:

1. Instead of a single Denisovan introgression into Australians 1353 generations ago of 4%, there are two introgressions 1353 generations ago into Papuans, one from D1 and one from D2

Hmm... Yes. Strictly speaking, the model described in Table S5 and Fig. S5 (fully described in S10.i, pg. 30--32) was implemented with a single time for the introgression of both Denisovan D1 and D2. They use this model to infer the admixture proportions of D1 and D2, assuming a total of 4% Denisovan ancestry. Just to be clear, I don't think the S5 model implements any parameters newly inferred by Jacobs et al., except the D1/D2 proportions.

In supplementary section S10.h (pg. 28--30) a separate analysis is done to determine the timing of the admixture. This doesn't use the demographic model, but instead uses the length distributions of denisovan-ancestry blocks assigned to D1 and D2. The admixture times I used are reported in the supplement at the top of page 30, but also in Fig 4, and the main text (page 6).

Do you think it's preferable to just stick to the Table/Figure S5 model (set T=1353*29=39237 for both introgressions), even though the one of the main results in the paper is that these introgressions occurred at different times?

[noscode]

That is very interesting question. I have not noticed this different times as I am trying to defeat test between qc model and yours :). I think that different times looks better.

I am not sure, may be I can just make changes in my qc PR? Anyway you have created another one.

[grahamgower]

I think its best if your QC doesn't touch the catalog model. If you're happy to use the two different introgression times, then feel free to submit your PR. If you rebase once #386 is merged, the tests should pass.

[noscode]

Okay I will wait. Test failure looks like a a bug actually.

[noscode]

Oh wait, I forgot about one more thing. You mixed up populations right after 1334 line:

- matrix_index=(pop["Papuan"], pop["Ghost"])),
+ matrix_index=(pop["CHB"], pop["Ghost"])),

- matrix_index=(pop["Ghost"], pop["Papuan"])),
+ matrix_index=(pop["Ghost"], pop["CHB"])),

May be you could add it to PR#386?

[jeromekelleher]

Just wanted to say here: it's fantastic to see the QC process working so well! Thanks @noscode and @grahamgower!

[noscode]

See PR#387.

[grahamgower]

Thanks @noscode. I owe you one QC.

This issue can be closed now.

[noscode]

That's sweet, thank you, @grahamgower!