"Diagnostic" plots and stats for hybrid models in msprime

[ChrystelleDelord]

Dear msprime team,

I am still new to msprime, and I hope I did not get any concept wrong before writing this!

In order to check-up that my simulations do not display anomalies in ancestry and LD patterns through time when using the hybrid DTWF/Hudson approach in msprime , I would like to record and/or output relevant information from the resulting tree sequence. I do not know if this is the best way to go but, so far, I was thinking of doing the following things:

• Get preliminary information about coalescence rates and proportion of uncoalesced lineages, using the DemographyDebugger utility for my models, as shown in the documentation,
• Simulate my models while sampling individual at several points in time (for example, each 2 generations),
• From the resulting tree sequence, get the actual coalescence rate trajectory and the actual proportion of ancestral lineages through time (I am not sure how to do this),
• From each sampling time, compute the MAF distribution, with the Fst is >1 pops,
• From samples, compare the SFS (or JSFS) between sampling times with varying time gap (e.g. 2 and 4 generations ago, 4 and 6, 6 and 8, etc, and then 10 and 20, 20 and 30...) using the distribution of residuals or the euclidian distance, and follow these through time.

I suppose, this would also be good to compare all of this with a "full-DTWF" or a "full-Hudson" simulation of the same model, as in Nelson et al. (2021) or Bhaskar et al. (2015).
Finally, I was also planning of doing the same procedure with simulations combining a nonWF model from SLiM3, and pyslim (I had posted on the SLiM-discuss googlegroup and was given the idea of following the SFS through time to see how steadily it is evolving).

May I ask you, what you think of my strategy? Is it relevant?
If so, should I have a look at the tree sequence node tables to compute the actual coalescence rate and the proportion of ancestral lineages after the simulation?

Thank you so much in advance for your help,

Kind regards
Chrystelle

[jeromekelleher]

Coming up with a good set of "diagnostic plots" is a great question @ChrystelleDelord. I'm not quite sure personally what would be a good thing to include, but hopefully others can chime in here with some suggestions.

The lessons learned paper may provide some useful background here.

[ChrystelleDelord]

Dear Jerome,

Thank you very much! I had a look at the Ragsdale et al. paper and it is very interesting. I had a few questions about the scripts used there, especially the ones to compute Figure 2 and other values.

So, I saw that this script is used to compute, amongst others, the coalescence rates using the coal_rates() function -which seems to be a method or an attribute of an object of type demography.DemoGraph.

When looking at the code for demography.coalescence_rates() from the package demography, it seems like the rate is kind of an "expected" rate (lines 118-127) calculated from demographic parameters of the model (maybe I am wrong!)
Is there a way (and is it acceptable/accurate) to get an "observed" rate from the node tables of a tree sequence, using their "time" value?

Thanks a bunch!
Chrystelle

[jeromekelleher]

I think a lot of the functionality in those scripts is now either in msrpime's numerical utilities or the Deme's package

We are working on getting these empirical distributions at the moment in tskit: tskit-dev/tskit#2119

[ChrystelleDelord]

Thank you Jerome,

I will have a look at all this and come back here to explain how it went, in case someone else would be interested.
See you soon,

Chrystelle

[jeromekelleher]

Thanks Chrystelle, that would be very helpful