Individual admixture computation

[alxsimon]

Hi everyone,

I am new to treesequences and started playing with some simulations where I create the model in Demes -> run in Slim -> produce TreeSequences and recapitate.
For instance I did a simple simulation of a 2 populations split with divergence, followed by the formation of an admixed population in the last 100 generations.
For instance:

description: Admixture between two populations - Sim 1 no migration
time_units: generations
generation_time: 30
demes:
- name: ancestral
  epochs:
  - {end_time: 100000.0, start_size: 5000}
- name: Pop_A
  ancestors: [ancestral]
  epochs:
  - {end_time: 100, start_size: 100}
  - {end_time: 0, start_size: 10000}
- name: Pop_B
  ancestors: [ancestral]
  epochs:
  - {end_time: 100, start_size: 100}
  - {end_time: 0, start_size: 10000}
- name: Pop_C
  start_time: 100
  ancestors: [Pop_A, Pop_B]
  proportions: [0.4, 0.6]
  epochs:
  - {end_time: 0, start_size: 100}

One of the thing I have a hard time to figure out right now is how to compute the admixture composition of each individual.
Could you point me to resources or tutorials that exist on that?
I feel like the introgression tutorial is close to what I want but I did not find an option in Slim equivalent to the msprime record_migration.

Or let me know if this is not possible at this time in tskit without using mutations to identify parental populations.

Thanks a lot,
Alexis

[petrelharp]

Generally, I think the thing to do is (a) Remember everyone in your tree sequence at the time you want to define "ancestral populations" for the purposes of visualizing admixture; and then (b) find IBD segments from your final sample back to those ancestors. Is this what you'd say, @gtsambos? (The documentation for this is not fleshed out yet, apologies.) This will take longer if the reference time is longer ago (more segements to track).

There ought to be a way to just output a summary instead of all the segments, but I don't see how to do it with our current set of statistics. It'd be fun to work out a way to do that.

[alxsimon]

Thanks a lot Peter!
I think I understand (hopefully) and this is what I ended up with:

def get_admixture_proportions(ts, admix_ind, popA_nodes, popB_nodes):
	max_time = ts.max_root_time - 1
	proportions_A = []
	for i in admix_ind:
		span_A = ts.ibd_segments(
			between=[popA_nodes, ts.individual(i).nodes],
			max_time=max_time,
		).total_span
		span_B = ts.ibd_segments(
			between=[popB_nodes, ts.individual(i).nodes],
			max_time=max_time,
		).total_span
		proportions_A.append(span_A/(span_A + span_B))
	return np.array(proportions_A)

popA_nodes = np.array([ts.individual(i).nodes for i in ts.individuals_alive_at(100, population=1)]).flatten()
popB_nodes = np.array([ts.individual(i).nodes for i in ts.individuals_alive_at(100, population=2)]).flatten()
admix_ind = ts.individuals_alive_at(0, population=3)

admixture_proportions = get_admixture_proportions(ts, admix_ind, popA_nodes, popB_nodes)

This seems to give me what I'm looking for, but I'll need to check in more detail if that's right.

For now it should do the trick as I don't consider a lot of individuals, but it would indeed be great if at some point there is a function to do that quickly.

[petrelharp]

Nice, yes, this is what I had in mind! I've not run it myself, but the only thing I notice is that the max_time can be the time of the reference individuals (so, 100 in your example). Or, a little older if you've a nonWF model, to make sure you get the birth times of all those. That should help with runtime (maybe a lot).

... oh wait, hold on. Looking at the documentation for ibd_segments, you've got to be a little careful here. Since everyone is IBD everywhere, I suspect that your method will actually return that everyone is 100% of both populations at the same time. Setting max_time to be equal to the reference time I think will fix this (since the only way two samples can be IBD at a time no longer ago than the birth time of one is if one is the descendant of the other). However, I was actually thinking of the link_ancestors function, which does what you want to do more directly here. But this is still just a method of the Table Collection - @gtsambos, what's the status of that?

[alxsimon]

Ok, not sure I grasp everything here.
That's what I understood from the ibd_segments, but using max_time 100 gives me 0 segments, in contrary to generation 101.
Starting from this, the span is increasing until the max_root_time which correspond indeed to everyone being 100% of both populations.

See for instance here where x is the max_time and y is the span of ibd for one sampled admixed individual (blue to A population, and orange to B population):

Regarding the link_ancestors function, I get an empty table when using one of the ancestral pop as ancestors and the admixed individuals as samples.

Anyway, I think this MWE is equivalent to my Slim simulation https://gist.github.com/alxsimon/32e6016153c50bc861fda98f4ae23c56

[petrelharp]

using max_time 100 gives me 0 segments, in contrary to generation 101.

That is not surprising to me. 101 must be the right number.

Those plots make sense! So, I think what you want is what you get from using max_time=101.

As for link_ancestors - could you post code that you used?

[alxsimon]

Yes that should work.
I need to recheck that but I think I had better correlation between allele freq based admixture estimation and ibd method when using the highest max_time possible than when using the 101 time, I think this is what was bothering me to start with.

For the link_ancestors you can have a look at the gist above https://gist.github.com/alxsimon/32e6016153c50bc861fda98f4ae23c56, but basically I used

popA_nodes = ts.samples(population=1)
popB_nodes = ts.samples(population=2)
admix_samples = ts.samples(population=3)

print(ts.tables.link_ancestors(admix_samples, popA_nodes))
print(ts.tables.link_ancestors(admix_samples, popB_nodes))

So that pops A and B are the samples at time 100.

[petrelharp]

Oh! I see what's going on. For link_ancestors to work, or for ibd with max_time = 101 to do what we want, we want to have nodes in all ancestral lineages at the reference time. That's what I mean by a "reference time" - things like admixture are not well-defined, and you have to give a particular time that it's defined with respect to. So, you therefore need to know, for each chunk of genome, who inherited from at that particular time. Say that's 101 generations ago. In SLiM you do that by doing RememberIndividuals on everyone alive at the time, so that any lineage from anyone later in time has to pass through one of those individuals. In msprime, to do this you need to add a census event.

What you're doing here is not sampling everyone in the population at the reference time, but instead taking a random sample of the individuals at that time. Since coalescence simulations are effectively infinite, these do not correspond to the lineages of your sample; but rather, they take some time to coalesce with them. From your plots, it looks like they coalesce over the next 100 generations or so. One could also define something analogous to admixture in this way (e.g., "on what proportion of the genome does the genome coalesce with one of these ancient samples from A before an ancient sample from B") but that's harder to think about and harder to compute.

[alxsimon]

Oooh I see! got it now.
This works great with the census event in the msprime simulation and I have the same result both with link_ancestors and the ibd_segments methods.
For the Slim simulation, I did RememberIndividuals in my previous simulations but the trick is that I used it in the late part of the generation when the admixed population is created. Switching it to the early part did the trick!

Thanks a lot for helping me understand what was going on, I'll try to write up a function working in my case and post it here as an answer.

[alxsimon]

Ok so thanks to @petrelharp, there are two ways of doing this.
The important point is to remember all ancestors at the point we define reference populations in our simulation, either using treeSeqRememberIndividuals in SLiM (I needed to use it in the early event in my case), or using a census event in msprime (not a sampling event).

Either find the admixture compositions with the ibd_segments (much slower):

def get_admixture_proportions(ts, admix_inds, popA_nodes, popB_nodes, max_time):
	proportions_A = []
	for i in admix_inds:
		span_A = ts.ibd_segments(
			between=[popA_nodes, ts.individual(i).nodes],
			max_time=max_time,
		).total_span
		span_B = ts.ibd_segments(
			between=[popB_nodes, ts.individual(i).nodes],
			max_time=max_time,
		).total_span
		proportions_A.append(span_A/(span_A + span_B))
	return np.array(proportions_A)

Or with the link_ancestors function (much faster):

def get_admixture_proportions_2(ts, admix_inds, popA_nodes, popB_nodes, max_time):
	proportions_A = []
	admix_nodes = np.array([ts.individual(ind).nodes for ind in admix_inds]).flatten()
	edges = ts.tables.link_ancestors(admix_nodes, np.concatenate((popA_nodes, popB_nodes)))
	for ind in admix_inds:
		child_nodes = ts.individual(ind).nodes
		A_edges = edges[np.isin(edges.child, child_nodes) & np.isin(edges.parent, popA_nodes)]
		B_edges = edges[np.isin(edges.child, child_nodes) & np.isin(edges.parent, popB_nodes)]
		span_A = np.sum(A_edges.right - A_edges.left)
		span_B = np.sum(B_edges.right - B_edges.left)
		proportions_A.append(span_A/(span_A + span_B))
	return np.array(proportions_A)

This gives exactly the same results and is equivalent to using fixed mutations to identify ancestral populations in my use case.

The ibd_segments version might be made faster if you get all the segments in one go and remember them all to filter them afterwards like in the edgeTable case. But I hadn't time to play with that.