Polish dog/wolf demography from #905

[gregorgorjanc]

I have taken the code from the draft PR at #905 and:

• removed all Boxer related code (as agreed in the draft PR add dog demographic model #905 (comment))
• added mutation rate 1e-8 in the model (which differs from the defined stdpopsim dog species rate of 4e-9)
• polished docs

I volunteer to QC this model to learn the process.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 99.84%. Comparing base (d6e0406) to head (6c49f26).
Report is 7 commits behind head on main.

Additional details and impacted files

@@           Coverage Diff           @@
##             main    #1632   +/-   ##
=======================================
  Coverage   99.84%   99.84%           
=======================================
  Files         130      132    +2     
  Lines        4416     4496   +80     
  Branches      608      608           
=======================================
+ Hits         4409     4489   +80     
  Misses          3        3           
  Partials        4        4           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[gregorgorjanc]

@agladstein could you have a look if you agree with my changes? I think they are straightforward. I have also been doing QC offline against the publication.

@nspope or @petrelharp could you review and hopefully merge - I have QC code ready and will share that in a separate PR.

[petrelharp]

Hm - we're missing something to get this model hit by tests...

[agladstein]

LGTM. I notice a few difference between your version and my last version - but, I think you are correct in your updates. (Looks like I had the wolf pop sizes swapped).

[gregorgorjanc]

Hm - we're missing something to get this model hit by tests...

I think I found the reason! Let's see what GitHub actions say.

[gregorgorjanc]

@agladstein I get this error now that tests run

E msprime._msprime.LibraryError: Attempt to move a lineage into an inactive population

I suspect this is because of how we removed Boxer - this line

model.add_population_split(time=T_ancDOG1, derived=["BOX", "BSJ"], ancestral="ancDOG1")

became

model.add_population_split(time=T_ancDOG1, derived=["BSJ"], ancestral="ancDOG1")

Do you have a suggestion on how to handle this? Perhaps we remove ancDOG1 and change BSJ Ne using model.add_population_parameters_change() or some other approach?

[agladstein]

Yes - that make sense. I was considering pointing that out, but I thought based on the msprime docs that it would still work.
I think your suggestion of removing ancDOG1 and change BSJ Ne using model.add_population_parameters_change() makes sense (and was what I was going to suggest before I checked the docs and thought keeping the split would still work).

If that doesn't work, I can try playing with it locally (but, I'm on a new computer and don't have my environments properly set up yet).

[gregorgorjanc]

To address this issue, I have now used the trunk model example suggestion from https://tskit.dev/msprime/docs/stable/demography.html#trunk-population-models

This now works on my end

import stdpopsim
species = stdpopsim.get_species("CanFam")
for x in species.demographic_models:
    print(x.id)
# EarlyWolfAdmixture_6F14

model = species.get_demographic_model("EarlyWolfAdmixture_6F14")
print(model.model)
# Demography(
# populations=[
# Population(initial_size=2639, growth_rate=0, name='BSJ', description='Basenji', extra_metadata={}, default_sampling_time=None, initially_active=True, id=0),
# Population(initial_size=1914, growth_rate=0, name='DNG', description='Dingo', extra_metadata={}, default_sampling_time=None, initially_active=None, id=1),
# Population(initial_size=5426, growth_rate=0, name='CHW', description='Chinese wolf', extra_metadata={}, default_sampling_time=None, initially_active=None, id=2),
# Population(initial_size=26092, growth_rate=0, name='ISW', description='Israeli wolf', extra_metadata={}, default_sampling_time=None, initially_active=None, id=3),
# Population(initial_size=11427, growth_rate=0, name='CRW', description='Croatian wolf', extra_metadata={}, default_sampling_time=None, initially_active=None, id=4),
# Population(initial_size=19446, growth_rate=0, name='GLJ', description='Golden jackal', extra_metadata={}, default_sampling_time=None, initially_active=None, id=5),
# Population(initial_size=1999, growth_rate=0, name='ancDOG', description='Ancestral ((Boxer, Basenji), Dingo)', extra_metadata={}, default_sampling_time=12795, initially_active=False, id=6),
# Population(initial_size=1393, growth_rate=0, name='ancWLF1', description='Ancestral (Israeli, Croatian) wolf', extra_metadata={}, default_sampling_time=13389, initially_active=False, id=7),
# Population(initial_size=12627, growth_rate=0, name='ancWLF', description='Ancestral ((Israeli, Croatian), Chinese) wolf', extra_metadata={}, default_sampling_time=13455, initially_active=False, id=8),
# Population(initial_size=44993, growth_rate=0, name='ancDW', description='Ancestral (dog, wolf)', extra_metadata={}, default_sampling_time=14874, initially_active=False, id=9),
# Population(initial_size=18169, growth_rate=0, name='root', description='Ancestral ((dog, wolf), golden jackal)root', extra_metadata={}, default_sampling_time=398262, initially_active=False, id=10)],
# events=[
# PopulationParametersChange(time=12102, initial_size=793, growth_rate=None, population='BSJ'),
# PopulationSplit(time=12795, derived=['BSJ', 'DNG'], ancestral='ancDOG'),
# PopulationSplit(time=13389, derived=['ISW', 'CRW'], ancestral='ancWLF1'),
# PopulationSplit(time=13455, derived=['ancWLF1', 'CHW'], ancestral='ancWLF'),
# PopulationSplit(time=14874, derived=['ancDOG', 'ancWLF'], ancestral='ancDW'),
# PopulationSplit(time=398262, derived=['ancDW', 'GLJ'], ancestral='root')],
# migration_matrix=array([[0.  , 0.  , 0.  , 0.18, 0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.  ],
# ...

contig = species.get_contig("chr38", mutation_rate=model.mutation_rate)
samples = {"BSJ": 2, "DNG": 2}
engine = stdpopsim.get_engine("msprime")
ts = engine.simulate(model, contig, samples)
print(ts.num_sites)
# 36158

@agladstein please check if these changes make sense.

[gregorgorjanc]

Hmm, I see that tests on GitHub fail though!? See https://github.com/popsim-consortium/stdpopsim/actions/runs/12556941637/job/35008826399?pr=1632

FAILED tests/test_models.py::TestCanFamEarlyWolfAdmixture_6F14::test_simulation_runs - msprime._msprime.LibraryError: Attempt to move a lineage into an inactive population

Any idea what is causing this error? When I install locally from my git repo fork, I don't get this error as shown above at #1632 (comment)

[petrelharp]

Splitting two pops into one is better practice than having a 'trunk' lineage in this case, since AFAIK neither BOX nor BSJ are a natural stand-in for ancestral dog. So, that's not the issue.

[gregorgorjanc]

@petrelharp we can’t include BOX because the publication didn’t estimate BOX parameters:( So, I have to use this trunk trick, I think. Happy to do something else, suggestions welcome.

[agladstein]

Not sure why you don't get the error locally, but GitHub fails. I was hoping to play with it locally myself, but I'll be out for next 2 weeks.
Possibly someone has a more technical explanation of why GitHub is failing, but I'd recommend trying your original idea of removing ancDOG1 and change BSJ Ne using model.add_population_parameters_change(), so it's just one branch with a population size change.

[petrelharp]

Ah, right, sorry - now I've read the information above. And, after some digging, I know what's going on. It's a bit of a mess.

Here's the reason for the error: if I insert print(samples) into the definition of test_simulation_runs( ) then I get

{'BSJ': 2, 'DNG': 2, 'CHW': 2, 'ISW': 2, 'CRW': 2, 'GLJ': 2, 'ancDOG': 2, 'ancWLF1': 2, 'ancWLF': 2, 'ancDW': 2, 'root': 2}

Digging in more, the problem is the samples from GLJ, not any of the ancestral populations. And... ah-ha. That's because we're starting migration between GLJ and ancDW before ancDW actually exists. Fixed!

(gee, that was annoying to track down - I kept getting distracted by thinking it was to do with getting samples from ancestral pops)

[petrelharp]

Great, now the docs build is failing with this (reproduced locally):

Exception occurred:
  File "/home/peter/micromamba/envs/stdpopsim/lib/python3.12/site-packages/docutils/writers/_html_base.py", line 852, in visit_entry
    if node.parent.parent.parent.stubs[node.parent.column]:
       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^
IndexError: list index out of range
The full traceback has been saved in /tmp/sphinx-err-mybaprij.log, if you want to report the issue to the developers.
Please also report this if it was a user error, so that a better error message can be provided next time.
A bug report can be filed in the tracker at <https://github.com/sphinx-doc/sphinx/issues>. Thanks!

Here's the full log:
sphinx-err-mybaprij.log

I've not got the bandwidth to debug this at the moment; has anyone else?

[petrelharp]

Never mind - it's because the parameter tables are (unquoted) CSV but some of these parameters have commas in their names:

Population size,"793",Ancestral (Boxer, Basenji) pop. size

should be

Population size,"793","Ancestral (Boxer, Basenji) pop. size"

[petrelharp]

I'm going to merge this! Could someone open the QC issue, please?

[gregorgorjanc]

Thanks @petrelharp!