Parameter estimates seem pretty far off

[Ellie Faber]

Hi Ryan,

I am running a series of simple models (SI, IM, AM, & SC) on an autosomal VCF dataset. Although my residuals for the best fit model (SC) look decent (at least to me), I am getting really strange estimates of divergence time for what we know about this system. These species are estimated to have split between 3-4 million years ago based on precious phylogenetic analyses, so I am not sure if I am just calculating the L parameter incorrectly or if something deeper is wrong.

This is how I am calculating L:

L = genome size * (number of SNPs entering analysis / number of SNPs that COULD enter analysis)

From the denominator, I've excluded sites in the VCF that were previously masked due to not passing filters or being masked as repeats. This number is therefore just the number of SNPs in the VCF being used. The numerator is the number of segregating sites in the jSFS for the two species in the analysis. 

Although my runs are converging, I am getting a theta value of around 172,470, and am using a mutation rate of 2.54e-08 and a generation time of 3.3 years, which results in the following estimates of divergence times:

Nref = 172,470/(4*2.54e-08*80,612,907)

Nref = 21057.9

T1 = 2.3

T1 (scaled) = 2*3.3*21057.9*2.3 = 315489.6

T2 = 1.0002

T2 (scaled)= 2*3.3*21057.9*1.0002= 139010

Td = T1+T2 = 454499.6

This seems extremely low given our knowledge of this system. I've attached a plot of my residuals + results of the best fit model. Any insight would be greatly appreciated!

[Ryan Gutenkunst]

Hello Ellie,

I am concerned that your model fit isn’t actually that good. The very large >500 residuals are swamping the signal in the plot. But just comparing the model and the data you can see qualitatively that the model is predicting many more fixed differences in pyres at intermediate frequencies in stephensi than the data show.

Your time conversion looks correct to me.

A challenge is that 3-4 million years ago with a generation time of 3.3 years is roughly 10^6 generations. If Ne is really of the order 10^4, then that divergence time is roughly 100 in population-size scaled units. In general, pop gen analysis isn’t going to have power unless that is at most 1-10. Over 100 time units, either all shared polymorphism would be lost, or if there’s migration it will have equilibrated, so there isn’t really signal to go that deep in time.

So to have power to resolve a 3-4 Mya divergence, your Ne needs to be roughly 100x what you’ve estimated here. I’m not sure that’s reasonable for rattlesnakes, since Ne ~ 10^6 is like Drosophila.

Biologically, maybe there is a long-ago phylogenetic divergence, but strong gene flow recently or long continuous gene flow has kept the polymorphism shared?

Best,

Ryan

--
You received this message because you are subscribed to the Google Groups "dadi-user" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dadi-user+...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/dadi-user/b0bbe4f8-2715-4455-8231-f27312c708b4n%40googlegroups.com.
<pyrr.steph.demo.plot.pdf>

[Ellie Faber]

HI Ryan,

Thank you for the quick response! That makes sense because knowing what we do about these species there appears to be more recent gene flow between them at their contact zone, so I wonder if that is impacting the patterns of shared polymorphism I'm seeing. For a divergence time this deep, do you have any suggestions for an alternative approach that might allow for more reasonable inference?

[Ryan Gutenkunst]

Hello Ellie,

The alternative model I’d study would be one in which you have two populations with no shared variation between them, then gene flow begins sometime in the past. Unfortunately, dadi doesn’t natively support that sort of model. (Although it could be hacked together.) The easiest approach would be to force it somewhat by first running a two_pop integration for a long T (maybe 20) with zero migration. That will be computationally slow, but maybe acceptable for your usage.

It’s surprising the phylogenetic estimate would be so deep with such strong shared polymorphism between the species.

Best,

Ryan

To view this discussion visit https://groups.google.com/d/msgid/dadi-user/c78ff8cb-be08-4514-a471-b96ce34ebd45n%40googlegroups.com.