Stdpopsim: QC for BosTau/HolsteinFriesian_1M13

[iechevarriaz]

PR for new model: https://github.com/popsim-consortium/stdpopsim/pull/558
Original paper: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3748359/
Parameter values: {reference to where parameter values can be found}
Population indices: {description of population indices}

Potential issues that might lead to differences between the model implementations:
{description of any potential issues}

QC'er requests: {tag potential developers to QC the model}

[gregorgorjanc]

Just a note here that I have "eyeballed" with @iechevarriaz the effective population sizes at different times from this paper https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3748359/ Figure 4b. These are the best estimates, from my knowledge.

[iechevarriaz]

Good morning! There are values for Ne and generations ago in the following supplementary material https://academic.oup.com/mbe/article/30/9/2209/1002512#supplementary-data

[agladstein]

I will work on this QC!

[grahamgower]

Any updates here @agladstein? It would be great to get Bos taurus in. Should we split the species and the demographic model into two commits?

[grahamgower]

I've taken a look at the demographic model while making pr #600, and it looks like the model implementation is not piecwise-constant, as in MacLeod et al. It's probably better to stick with the publication here, unless there's a good reason not to?

I.e. at the moment, the model has:

        demographic_events=[
            msprime.PopulationParametersChange(
                time=1,
                initial_size=90,
                growth_rate=-0.095894024,
                population_id=0,
            ),  # Ne 90 to 120
            msprime.PopulationParametersChange(
                time=4, growth_rate=-0.24465639, population_id=0
            ),  # Ne 120 to 250
            msprime.PopulationParametersChange(
                time=7, growth_rate=-0.0560787, population_id=0
            ),  # Ne 250 to 350

But MacLeod et al. actually shows something more like:

        demographic_events=[
            msprime.PopulationParametersChange(
                time=1,
                initial_size=90,
                growth_rate=0,
                population_id=0,
            ),  # Ne = 90
            msprime.PopulationParametersChange(
                time=4, initial_size=120, population_id=0
            ),  # Ne = 120
            msprime.PopulationParametersChange(
                time=7, initial_size=250, population_id=0
            ),  # Ne = 250

[gregorgorjanc]

Thanks @grahamgower! @iechevarriaz and I were not sure what would be best. In reality the change was likely gradual, but I see merit in exactly following the publication. Please advise.

[grahamgower]

Well using a growth rate like that will make each epoch have an exponential population size change, which is probably not what is desired either. So I would recommend just sticking with a piecewise-constant model. (Maybe in the future we'll try to include demographic models with nice smooth changes through time as in https://www.biorxiv.org/content/10.1101/2020.09.21.307355v1)

[gregorgorjanc]

OK, makes sense. This feedback is awesome.

[agladstein]

@grahamgower just to clarify, we still need to do a blind implementation for QC of this model right?

[grahamgower]

Yes please @agladstein! This will also need #601 and #606 before QC is complete. Maybe you'd prefer to work from the values in the parameter table (#606) and double check they match the figure in the paper?

[gregorgorjanc]

Yes please @agladstein! This will also need #601 and #606 before QC is complete

@agladstein and @grahamgower #609 addresses the #601 and #606

[agladstein]

Great. I'll wait until https://github.com/popsim-consortium/stdpopsim/pull/609 gets merged.

[gregorgorjanc]

@agladstein #609 has been merged in. Maybe we do QC together at the "Grow the Zoo" hackathon?

[agladstein]

I was planning to add in a dog demographic model.
Maybe there is time to do both?

[gregorgorjanc]

Purfect! You can QC my cattle species&model and I QC for your dog species&model.