split vs. trunk model differences

[mshpak76]

For a model of three subpopulations (A, (B,C)), I am trying to understand the difference between modeling it with ancestral splits (i.e. explicit ancestor BC for (B,C) and ABC for (A,BC)) vs. the "trunk" model where C merges with B and B merges with A. In the absence of bottlenecks at the time of the splits, shouldn't this give identical results? At least that is what I thought until I tested the two models using the same mutation rate etc, i.e. split model, i.e. common ancestor ABC splits into A and BC, while BC splits into B,C

demography = msprime.Demography()
demography.add_population(name = "ABC", initial_size = 2e+6)
demography.add_population(name = "BC", initial_size = 2e+5)
demography.add_population(name = "A", initial_size = 2e+6)
demography.add_population(name = "B", initial_size = 4e+5)
demography.add_population(name = "C", initial_size = 4e+5)
demography.add_population_split(time = 4e+4, derived = ["B","C"], ancestral = "BC")
demography.add_population_split(time = 4e+5, derived = ["A", "BC"], ancestral = "ABC")

vs. trunk model where A is ancestral to B which is ancestral to C

demography = msprime.Demography()
demography.add_population(name = "A", initial_size = 2e+6)
demography.add_population(name = "B", initial_size = 4e+5)
demography.add_population(name = "C", initial_size = 4e+5)
demography.add_population_split(time = 4e+4, derived = "C", ancestral = "B")
demography.add_population_split(time = 4e+5, derived = "B", ancestral = "A")

For the second "trunk" model, I consistently get distance(A,B) < distance(A,C), whereas they are essentially identical for the first (split) model. This is particularly apparent when I calculate Fst, i.e. Fst(A,B) < Fst(A,C). It is not apparent to me why there should be any differences at all. As far as I can tell, the only difference between the models is whether I label the ancestral lineage of (B,C) as "B" or as BC, since any mutations accumulated along that lineage will be inherited by B and C regardless.

[grahamgower]

Hey @mshpak76. The two ways of specifying a model should give equivalent results. One problem with your comparison is that in the first model, the population size of BC is 2e5, but in your second model the population size of B is 4e5 (even after C merges into B).

[mshpak76]

That was a typo when changing parameters from two different scripts, BC should be 4e5 just like B (I've corrected it above). The latter model still consistently gives an asymmetry between B and C whereas the former does not.

I think that the problem is that I need to set both A and B to initially_active for the latter model to be consistent with the first (otherwise there is a gap in the time evolution of the B lineage in comparison to C)

[grahamgower]

Whether populations are initially_active shouldn't affect things---an inactive population can't be sampled from and events involving inactive populations will give an error (to help catch model misspecification errors).

Can you post your full script?

[mshpak76]

Here's an example of code to generate the (A,(B,C)) genealogy that gives Fst(A,C) > Fst(A,B) (on average), even though the two Fst values should be approximately equal. Strangely enough, if I add initially_active = True to the add_population for A and B, I get nearly equal Fst for the A,C and A,B pairs.

My question is why the demographic model below should result in any asymmetry of A wrt B and C:

samp_size = 25
demography = msprime.Demography()
demography.add_population(name = "A", initial_size = 20000)
demography.add_population(name = "B", initial_size = 4000)
demography.add_population(name = "C", initial_size = 4000) 
demography.add_population_split(time=400, derived=["C"], ancestral="B")
demography.add_population_split(time=4000, derived=["B"], ancestral="A")

def run_sim(num_samples, ancestry_seed, mutation_seed):
   ancestral_ts = msprime.sim_ancestry(samples={"A":samp_size, "B":samp_size, "C":samp_size}, demography= demography, sequence_length=10000, recombination_rate=1.14e-8, gene_conversion_rate=5.9e-8, gene_conversion_tract_length=1, random_seed= ancestry_seed)    
   mutated_ts = msprime.sim_mutations(ancestral_ts, rate=1.22e-8, random_seed=mutation_seed, model=msprime.BinaryMutationModel())
   return mutated_ts

mts = run_sim(num_replicates, *seeds[rep_index])

Note: There's an additional block of code that generates random seeds for the mutation and ancestry, as well as a loop to iterate mts and write output to a file.

[grahamgower]

Ah, I think what is happening is that B has initially_active=False because it's an ancestor in your split event, and this causes the default sampling time to be 400 (the time at which B becomes active). By setting initially_active=True, the default sampling time is 0. (And likewise for A)

[mshpak76]

That makes sense. That would explain Fst(A,B) < Fst(A,C), if C is sampled at 0 and B at 400.

…

On Fri, Jun 17, 2022 at 3:02 PM Graham Gower ***@***.***> wrote: Ah, I think what is happening is that B has initially_active=False because it's an ancestor in your split event, and this causes the default sampling time to be 400 (the time at which B becomes active). By setting initially_active=True, the default sampling time is 0. — Reply to this email directly, view it on GitHub <#2072 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AUZMZCNFHY5VV3SEPMYDKPDVPTKU3ANCNFSM5ZDFJTYA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

-- ======================= Max Shpak, Ph.D. Department of Genetics University of Wisconsin Madison, WI 53706