best way to model admixture in final generation

[mshpak76]

Rather than simulating gene flow in every generation, I would like to model a scenario where there was admixture in a single generation (perhaps in the generation before sampling), with no migration during the rest of the population's history. Is there a straightforward way to do this? Using

demography.set_migration_rate(source = "p1", dest = "p2", rate = mrate)

would just have gene flow in each generation, so I'm not sure how to simulate admixture in a specific generation using msprime. I did find a function add_admixture, but I'm not sure if the following gives me the desired result of admixture a generation prior to sampling (mrate percent of population p1 having ancestry in p2 in the generation prior to sampling).

add_admixture(time=1, derived = "p1", ancestral = "p2", proportion = mrate)

I'm also not sure what it means to set the ancestral population as active vs. inactive with this function. If I have a genealogy with 2 splits (generating populations p1, p2, p3) and implement admixture in previous t=1, I'm presumably still able to sample the genotypes in generation 0?

Additionally, I get the following error message in association with add_admixture, even though I use the same names/format of input arguments as for add_population_split (which executes without issues):

demography.add_population_split(time=2000, derived=["B"], ancestral="A")
demography.add_admixture(time=1, derived=["C"], ancestral="B", proportions = 0.2)

The error message that I get for add_admixture is

"_KeyError: "Keys must be either string population names or integer IDs:identifier '['C']' is of type <class 'list'>"_, 

which (again) doesn't occur when the same population names are used for add_population

[petrelharp]

I think you're working with this mental model:

----------------------------------------- A
                                      \
----------------------------------------- B
   -------> time ---->                 t

i.e., that "at time t population B got a bunch of ancestry from A". However, msprime is steering you towards this mental model:

----------------------------------------- A
                                      \
                                       --- C
                                      /
----------------------------------------- B
   -------> time ---->                 t

i.e., "at time t admixture between A and B created population C". In this picture, both A and B are still around. To do something equivalent to the first picture, you'd have B no longer be around, i.e., make B inactive at the time of admixture:

----------------------------------------- A
                                      \
                                       --- C
                                      /
--------------------------------------x B
   -------> time ---->                 t

To do something like the above, you'd want to not have a population C beforehand, then call

demography.add_admixture(time=1, derived=["C"], ancestral=["A", "B"], proportions = [0.2, 0.8])

So - your call to add_admixture doesn't make sense because there's only one ancestral population, and the proportions given don't add up to 1.

[mshpak76]

In my model, populations B,C split 1K generations ago while (B,C) split from A 2K generations ago. So at one generation prior to sampling, all three populations are present. I would like to simulate the first scenario in your schematic, specifically a case where 20% of individuals in generation 1 are migrants from B while the remaining 80% are C residents, at the generation prior to sampling. I tried to do this with:

demography.add_admixture(time=1, derived=["C"], ancestral=["B", "C"], proportions = [0.2, 0.8])

with C ancestral to the remaining 80% C but I got an error message stating that the admixture event was out of chronological order with respect to the other demography functions (i.e. the splits). Surely there is some way to model the admixture scenario that I'm trying to simulate, i.e. with an established ( (B,C), A) tree established and admixture of C from B in the generation prior to sampling (or indeed any specified generation between 0 and 1K, the time of the B,C split)?

You mentioned that to do this, I would need to make (in my example) A inactive. However, at time 0, I want to sample from A,B,C, and I don't want all of C to be derived from B, only 0.2 (and various other values < 1), so I don't think that add_admixture is going to work here, but I'm guessing that something else might? One rather contrived method might be to create a "shell" population D that splits from C in generation 1 from an 80%, 20% C,B admixture, but I'm hoping that there's a more elegant way to do this.

Thank you, Max

[petrelharp]

Okay, let's see. I'm guessing the error you got was

ValueError: Events must be time-sorted. Please use demography.sort_events() if you add events out of order.

This is fixed if you either (1) do as the error message suggests: run demography.sort_events(), or (2) rearrange your code.

But, I think this does what you want:

demography = msprime.Demography()
demography.add_population(name="A", initial_size=100)
demography.add_population(name="BC", initial_size=100)
demography.add_population(name="B", initial_size=100)
demography.add_population(name="C", initial_size=100)
demography.add_population(name="today", initial_size=100)
demography.add_admixture(1, derived="today", ancestral=["B", "C"], proportions=[0.2, 0.8])
demography.add_population_split(1000, derived=['B', 'C'], ancestral='BC')
demography.add_population_split(2000, derived=['BC'], ancestral='A')

You should check this. This produces:

>>> print(demography)
Demography
╟  Populations
║  ┌────────────────────────────────────────────────────────────────────────────────────────────────┐
║  │ id │name   │description  │initial_size  │ growth_rate │  default_sampling_time│extra_metadata  │
║  ├────────────────────────────────────────────────────────────────────────────────────────────────┤
║  │ 0  │A      │             │100.0         │      0      │                  2e+03│{}              │
║  │ 1  │BC     │             │100.0         │      0      │                  1e+03│{}              │
║  │ 2  │B      │             │100.0         │      0      │                      0│{}              │
║  │ 3  │C      │             │100.0         │      0      │                      0│{}              │
║  │ 4  │today  │             │100.0         │      0      │                      0│{}              │
║  └────────────────────────────────────────────────────────────────────────────────────────────────┘
╟  Migration Matrix
║  ┌────────────────────────────────┐
║  │       │ A │ BC │ B │ C │ today │
║  ├────────────────────────────────┤
║  │      A│ 0 │ 0  │ 0 │ 0 │   0   │
║  │     BC│ 0 │ 0  │ 0 │ 0 │   0   │
║  │      B│ 0 │ 0  │ 0 │ 0 │   0   │
║  │      C│ 0 │ 0  │ 0 │ 0 │   0   │
║  │  today│ 0 │ 0  │ 0 │ 0 │   0   │
║  └────────────────────────────────┘
╟  Events
║  ┌─────────────────────────────────────────────────────────────────────────────────┐
║  │  time│type        │parameters          │effect                                  │
║  ├─────────────────────────────────────────────────────────────────────────────────┤
║  │     1│Admixture   │derived=today       │Moves all lineages from admixed         │
║  │      │            │ancestral=[B, C]    │population 'today' to ancestral         │
║  │      │            │proportions=[0.20,  │populations. Lineages move to 'B' with  │
║  │      │            │0.80]               │proba 0.2; 'C' with proba 0.8. Set      │
║  │      │            │                    │'today' to inactive, and all migration  │
║  │      │            │                    │rates to and from 'today' to zero.      │
║  │┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈│
║  │  1000│Population  │derived=[B, C],     │Moves all lineages from derived         │
║  │      │Split       │ancestral=BC        │populations 'B' and 'C' to the          │
║  │      │            │                    │ancestral 'BC' population. Also set     │
║  │      │            │                    │the derived populations to inactive,    │
║  │      │            │                    │and all migration rates to and from     │
║  │      │            │                    │the derived populations to zero.        │
║  │┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈┈│
║  │  2000│Population  │derived=[BC],       │Moves all lineages from the 'BC'        │
║  │      │Split       │ancestral=A         │derived population to the ancestral     │
║  │      │            │                    │'A' population. Also set 'BC' to        │
║  │      │            │                    │inactive, and all migration rates to    │
║  │      │            │                    │and from the derived population to      │
║  │      │            │                    │zero.                                   │

and using demesdraw like demesdraw.tubes(demography.to_demes()) gets

If you have any further issues, please post code, and also error messages.

[mshpak76]

Thanks, I used something close to what you suggested (other than using a single label B for the ancestral BC lineage) and it seems to work, insofar as not throwing any errors.

However, it's unclear to me how to now sample genotypes from the tree, as there are now 4 lineages "today", "C", "B", "A" rather than the original 3. Since I'm using "today" as a surrogate for C after admixture, I want to sample "today", "B", "A". However, the arrays seem to only consist of 3 populations rather than 4 at the time of sampling:

Demographic model:

demography = msprime.Demography()
demography.add_population(name = "A", initial_size = 20000, initially_active=True)
demography.add_population(name = "B", initial_size = 10000, initially_active=True)
demography.add_population(name = "C", initial_size = 10000)
demography.add_population(name = "today", initial_size = 10000)
demography.add_admixture(1, derived="today", ancestral=["B","C"], proportions=[0.2,0.8])
demography.add_population_parameters_change(time=1000, population="B", initial_size=5000)
demography.add_population_split(time=1000, derived=["C"], ancestral="B")
demography.add_population_split(time=2000, derived=["B"], ancestral="A")

Running the coalescent and sampling (this was originally in a replication loop for the code, to simplify I set the loop indices to 1):

num_replicates = 1
rep_index = 1
samp_size = 25

mts = run_sim(num_replicates, *seeds[rep_index])
seg_sites = mts.segregating_sites(sample_sets=[mts.samples(population=0), mts.samples(population=1), mts.samples(population=2),]).tolist()
diversity = mts.diversity(sample_sets=[mts.samples(population=0), mts.samples(population=1), mts.samples(population=2),]).tolist()
hap_size = 2*samp_size
Genotype_Array = mts.genotype_matrix()
print(len(Genotype_Array))
print(len(Genotype_Array[1,:]))

The second print statement gives me a Genotype_Array length of 150 (50 haplotypes x 3) rather than 200 (50 haplotypes x 4). I'm not sure what these 3 populations correspond to - are they "A","B","C" vs. "A","B", "today", or something else? Since "today" was generated from both B and C, I'm not sure which one I need to sample if I want the equivalent of post-admixture C with A,B as outgroups. Similarly, in my seg_sites and nucleotide diversity calculations, what do populations 0,1,2 now correspond to? Are they A,B,C or A,B,"today" (there is no 4th population , i.e. poulation=3 to sample).

[petrelharp]

Hm - I can't tell why you're getting 150 haplotypes because that's wrapped up in your run_sim code, which isn't here. (and instead of posting, better to simplify it down to just a call to msprime.sim_ancestry( ).)

To see which population is which; do

for p in ts.populations():
  print(p)

and there are many ways to do it, but to see which sample is from which population, you can do:

for n in ts.samples():
  print(ts.node(n))

Since "today" was generated from both B and C, I'm not sure which one I need to sample if I want the equivalent of post-admixture C with A,B as outgroups.

Why not C, since that's the one you want?

[mshpak76]

When I ran

for p in mts.populations():
     print(p)

I get the output

Population(id=0, metadata={'description': '', 'name': 'A'})
Population(id=1, metadata={'description': '', 'name': 'B'})
Population(id=2, metadata={'description': '', 'name': 'C'})
Population(id=3, metadata={'description': '', 'name': 'today'})

Suggesting that there are indeed 4 populations. However, if I substitute 3 for population 2 in:

seg_sites = mts.segregating_sites(sample_sets=[mts.samples(population=0), mts.samples(population=1), mts.samples(population=3),]).tolist()

I get an error message suggesting that population 3 ("today") doesn't exist or is empty:

ile "/opt/miniconda3/lib/python3.8/site-packages/tskit/trees.py", line 6363, in __one_way_sample_set_stat
raise ValueError("Sample sets must contain at least one element")
 ValueError: Sample sets must contain at least one element

The segregating_sites function works fine for populations 0,1,2 but not 3.

"Today" is the population I wish to sample because I'm using it as a surrogate for what C would be if in generation 1 20% of it had ancestry in population B. This would be much simpler if I could just model migration from B to C directly in a single generation, but if using this placeholder as a placeholder is the only way to do this, I need to figure out how to sample from it and why it's (apparently) empty.

[petrelharp]

Well, maybe you didn't take any samples from that population - what samples were specified? (This is hidden in the call to run_sim( ).)