Modeling population growth

[Madee Chase]

Dear Ryan,

 

I am trying to infer the demographic history of two hybridizing plants, and I had a couple questions about the models I’ve been fitting. We’ve fitted and compared several models already, including IM, SC, PSC, AM and SI. Additionally we’ve incorporated other parameters into each of these model families, such as heterogeneous migration rates and population growth.

 

I want to make sure we’re implementing these models correctly and that our interpretation of the time parameters are correct. Would you mind taking a look at the model for secondary contact, both with and without a period of population expansion?

 

Here is what I have for secondary contact without population expansion:

 

    # phi for the equilibrium ancestral population

    phi = dadi.PhiManip.phi_1D(xx)

    # Now do the divergence event

    phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)

    # We set the population sizes after the split to nu1 and nu2 and the migration rate to zero

    phi = dadi.Integration.two_pops(phi, xx, Ts, nu1, nu2, m12=0, m21=0)

    # We keep the population sizes after the split to nu1 and nu2 and set the migration rates to m12 and m21

    phi = dadi.Integration.two_pops(phi, xx, Tsc, nu1, nu2, m12=m12, m21=m21)

 

And then this is my model with expansion:

 

    # phi for the equilibrium ancestral population

    phi = dadi.PhiManip.phi_1D(xx)

    # Now do the divergence event

    phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)

    # Population growth

    nu1_func = lambda t: numpy.exp(numpy.log(nu1) * t/Tg)

    nu2_func = lambda t: numpy.exp(numpy.log(nu2) * t/Tg)

    # We set the population sizes after the split to nu1 and nu2 and the migration rate to zero

    phi = dadi.Integration.two_pops(phi, xx, Ts, nu1_func, nu2_func, m12=0, m21=0)

    # We keep the population sizes after the split to nu1 and nu2 and set the migration rates to m12 and m21

    phi = dadi.Integration.two_pops(phi, xx, Tsc, nu1_func, nu2_func, m12=m12, m21=m21)

 

 

In this case Ts = Scaled time between split and population growth

Tg = Scaled time between expansion and secondary contact

Tsc = Scaled time between secondary contact and the present

 

Is this correct?

 

I’ve also attached our spectrum. Intuitively does it give you a sense of what kind of history of gene flow you might expect?

[Gutenkunst, Ryan N - (rgutenk)]

Hello Madee,

See below.

On Nov 14, 2016, at 4:23 PM, Madee Chase <madee...@gmail.com> wrote:

Dear Ryan,

 

I am trying to infer the demographic history of two hybridizing plants, and I had a couple questions about the models I’ve been fitting. We’ve fitted and compared several models already, including IM, SC, PSC, AM and SI. Additionally we’ve incorporated other parameters into each of these model families, such as heterogeneous migration rates and population growth.

 

I want to make sure we’re implementing these models correctly and that our interpretation of the time parameters are correct. Would you mind taking a look at the model for secondary contact, both with and without a period of population expansion?

 

Here is what I have for secondary contact without population expansion:

 

    # phi for the equilibrium ancestral population

    phi = dadi.PhiManip.phi_1D(xx)

    # Now do the divergence event

    phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)

    # We set the population sizes after the split to nu1 and nu2 and the migration rate to zero

    phi = dadi.Integration.two_pops(phi, xx, Ts, nu1, nu2, m12=0, m21=0)

    # We keep the population sizes after the split to nu1 and nu2 and set the migration rates to m12 and m21

    phi = dadi.Integration.two_pops(phi, xx, Tsc, nu1, nu2, m12=m12, m21=m21)

Looks correct.

 And then this is my model with expansion:

 

    # phi for the equilibrium ancestral population

    phi = dadi.PhiManip.phi_1D(xx)

    # Now do the divergence event

    phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)

    # Population growth

    nu1_func = lambda t: numpy.exp(numpy.log(nu1) * t/Tg)

    nu2_func = lambda t: numpy.exp(numpy.log(nu2) * t/Tg)

    # We set the population sizes after the split to nu1 and nu2 and the migration rate to zero

    phi = dadi.Integration.two_pops(phi, xx, Ts, nu1_func, nu2_func, m12=0, m21=0)

    # We keep the population sizes after the split to nu1 and nu2 and set the migration rates to m12 and m21

    phi = dadi.Integration.two_pops(phi, xx, Tsc, nu1_func, nu2_func, m12=m12, m21=m21)

  

In this case Ts = Scaled time between split and population growth

Tg = Scaled time between expansion and secondary contact

Tsc = Scaled time between secondary contact and the present

I’m not sure this is quite what you’re looking for. I assume you want to set Tg = Ts + Tsc. That will model populations that split and start growing, then have secondary contact after Ts time units. Tg is set to Ts+Tsc so that the growth is scaled to match the total timespan after the split.

Is this correct?

 

I’ve also attached our spectrum. Intuitively does it give you a sense of what kind of history of gene flow you might expect?

It’s hard to have much intuition about gene flow just visually. It is interesting that you have a lot of SNPs that are high frequency in PUN but almost absent in AUS. That may be challenging to model.

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 post to this group, send email to dadi...@googlegroups.com.
Visit this group at https://groups.google.com/group/dadi-user.
For more options, visit https://groups.google.com/d/optout.
<ry_noDLR_85.pdf>

--

Ryan Gutenkunst
Assistant Professor of Molecular and Cellular Biology, University of Arizona
phone: (520) 626-0569, office: LSS 325, web: http://gutengroup.mcb.arizona.edu

Latest papers: 
“Selection on network dynamics drives differential rates of protein domain evolution”
PLoS Genetics; http://dx.doi.org/10.1371/journal.pgen.1006132
"Triallelic population genomics for inferring correlated fitness effects of same site nonsynonymous mutations"
Genetics; http://dx.doi.org/10.1534/genetics.115.184812
"Whole genome sequence analyses of Western Central African Pygmy hunter-gatherers reveal a complex demographic history and identify candidate genes under positive natural selection"
Genome Research; http://dx.doi.org/10.1101/gr.192971.115

[Madee Chase]

Thanks for the reply! Just to be clear, I would modify the growth function to this, correct?

nu1_func = lambda t: numpy.exp(numpy.log(nu1) * t/(Ts+Tsc))

[Gutenkunst, Ryan N - (rgutenk)]

Exactly.

--
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 post to this group, send email to dadi...@googlegroups.com.
Visit this group at https://groups.google.com/group/dadi-user.
For more options, visit https://groups.google.com/d/optout.

[Kun Wang]

Hi, Ryan,

I'm confused about the nu_func with multi stages.

Assume nu1_0 is the population size in the time of split, and nu1 is the population size in the present. 

Assume there are two stages, stage 1 with no gene flow (T1), stage 2 with gene flow (T2):

I want to know which one is right:

Model1:

T = T1+T2

nu1_func = lambda t: nu1_0 * (nu1/nu1_0)**(t/T1+T2)

nu2_func = lambda t: nu2_0 * (nu2/nu2_0)**(t/T1+T2)

phi = dadi.Integration.two_pops(phi, xx, T1, nu1_func, nu2_func, m12=0, m21=0)

phi = dadi.Integration.two_pops(phi, xx, T2, nu1_func, nu2_func, m12=m12, m21=m21)

Model 2:

T = T1+T2

nu1_func = lambda t: nu1_0 * (nu1/nu1_0)**(t/T)

nu2_func = lambda t: nu2_0 * (nu2/nu2_0)**(t/T)

phi = dadi.Integration.two_pops(phi, xx, T1, nu1_func, nu2_func, m12=0, m21=0)

nu1_func = lambda t: nu1_0 * (nu1/nu1_0)**((t+T1)/T)

nu2_func = lambda t: nu2_0 * (nu2/nu2_0)**((t+T1)/T)

phi = dadi.Integration.two_pops(phi, xx, T2, nu1_func, nu2_func, m12=m12, m21=m21)

I have tried to compare the other two model (Model3, Model4) with simple expansion model (Model 5):

Model 3:

T = T1+T2

nu1_func = lambda t: nu1_0 * (nu1/nu1_0)**(t/T1+T2)

nu2_func = lambda t: nu2_0 * (nu2/nu2_0)**(t/T1+T2)

phi = dadi.Integration.two_pops(phi, xx, T1, nu1_func, nu2_func, m12=m12, m21=m21)

phi = dadi.Integration.two_pops(phi, xx, T2, nu1_func, nu2_func, m12=m12, m21=m21)

Model 4:

T = T1+T2

nu1_func = lambda t: nu1_0 * (nu1/nu1_0)**(t/T)

nu2_func = lambda t: nu2_0 * (nu2/nu2_0)**(t/T)

phi = dadi.Integration.two_pops(phi, xx, T1, nu1_func, nu2_func, m12=m12, m21=m21)

nu1_func = lambda t: nu1_0 * (nu1/nu1_0)**((t+T1)/T)

nu2_func = lambda t: nu2_0 * (nu2/nu2_0)**((t+T1)/T)

phi = dadi.Integration.two_pops(phi, xx, T2, nu1_func, nu2_func, m12=m12, m21=m21)

I believe one of Model 3 and Model 4 will be equal to simple expansion model (Model 5). However, I found the performance is very different between Model 3/4 and Model 5.

Model 5:

T = T1+T2

nu1_func = lambda t: nu1_0 * (nu1/nu1_0)**(t/T)

nu2_func = lambda t: nu2_0 * (nu2/nu2_0)**(t/T)

phi = dadi.Integration.two_pops(phi, xx, T, nu1_func, nu2_func, m12=m12, m21=m21)

result of Model 3:

result of Model 4:

result of Model 5:

I want to know what should i do with expansion model of multiple stages.

Best wishes!

Yours,

Kun

[Gutenkunst, Ryan N - (rgutenk)]

Hello Kun,

Neither of your proposals is correct. One way to do it is:

T = T1+T2

nu1_func = lambda t: nu1_0 * (nu1/nu1_0)**(t/T)

nu2_func = lambda t: nu2_0 * (nu2/nu2_0)**(t/T)

phi = dadi.Integration.two_pops(phi, xx, T1, nu1_func, nu2_func, m12=0, m21=0)

nu1_temp = nu1_func(T1)

nu2_temp = nu2_func(T1)

nu1_func = lambda t: nu1_temp*(nu1/nu1_temp)**(T/T2)

nu2_func = lambda t: nu2_temp*(nu2/nu2_temp)**(T/T2)

phi = dadi.Integration.two_pops(phi, xx, T2, nu1_func, nu2_func, m12=m12, m21=m21)

Or you could do:

T = T1+T2

nu1_func = lambda t: nu1_0 * (nu1/nu1_0)**(t/T)

nu2_func = lambda t: nu2_0 * (nu2/nu2_0)**(t/T)

phi = dadi.Integration.two_pops(phi, xx, T1, nu1_func, nu2_func, m12=0, m21=0)

phi = dadi.Integration.two_pops(phi, xx, T2, nu1_func, nu2_func, m12=m12, m21=m21, initial_t=T1)

Best,

Ryan

[Kun Wang]

Hi Ryan,

Thank you for your advice.

However, I still got different figure when I use the these two new functions. The scripts and spectrum file is attached. Is there something wrong in my scripts?

Best wishes!

Kun

[Gutenkunst, Ryan N - (rgutenk)]

Hello Kun,

Two errors:

In model1: The second set of functions should be:

    nu1_func = lambda t: nu1_temp*(nu1/nu1_temp)**(t/T2)

    nu2_func = lambda t: nu2_temp*(nu2/nu2_temp)**(t/T2)

 (You had **(T/T2) instead)

In model2: The second integration should be:

    phi = dadi.Integration.two_pops(phi, xx, T1+T2, nu1_func, nu2_func,

                                    initial_t=T1)

 (My mistake)

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 post to this group, send email to dadi...@googlegroups.com.
Visit this group at https://groups.google.com/group/dadi-user.
For more options, visit https://groups.google.com/d/optout.

<dadi.plot.gz>

[Kun Wang]

Hello Ryan,

It works! Thank you very much!

Best wishes!

Kun