Identify segments from a given population after admixture

[leovallini]

Hi msprime community,
I'm sorry for the noob question but I was not able to find the solution elsewhere:

I have a simple population demography where population D is the admixture of population B and C with proportion 0.75 and 0.25 (see figure).

What I want to do is simulate multiple individuals from the admixed population D and for each of them(crucial, I want to know that for single individuals, not population-wise) save the left and right coordinates of the chunks that come from population C.
So far I only managed to do that for a single individual using record_migrations = True and selecting migration events where the source population is C, when I try to do that for more than one individual in the admixed D population, however, I get many migration events where the node does not coincide with one of my samples nodes; I think the reason for that is that the segments referred to in these migrations are shared by more than one individual in the admixed population, but how do I know which ones?

Below is the code for what I have done so far, which works fine if I simulate 1 individual from population D, but not when there is more than 1.

import msprime
import numpy as np

# define demography #
demography = msprime.Demography()
demography.add_population(name = "A", initial_size = 10000)
demography.add_population(name = "B", initial_size = 3000)
demography.add_population(name = "C", initial_size = 3000)
demography.add_population(name = "D", initial_size = 5000)
demography.add_admixture(1000, derived = "D", ancestral = ["B", "C"], proportions = [0.75, 0.25],)
demography.add_population_split(time = 2000, derived = ["B", "C"], ancestral = "A")

# simulate #
def run_simulation(sequence_length, demography):
    ts = msprime.sim_ancestry(
        demography = demography,
        sequence_length = sequence_length,
        samples = {"D": 10, "B": 10, "C": 10},
        ploidy = 1,
        recombination_rate = 1e-8,
        record_migrations = True
    )
    return ts

# identify all tracts that migrated backwards in time from D to C # 
# inspired by https://tskit-dev.github.io/tutorials/introgression.html #
def get_migrating_tracts(ts):
    C_id = [p.id for p in ts.populations() if p.metadata['name'] == 'C'][0]
    migrating_tracts = [("node", "left", "right")]
    for migration in ts.migrations():
        if migration.dest == C_id:
            migrating_tracts.append((migration.node, migration.left, migration.right))
    return np.array(migrating_tracts)

###################################################################

ts = run_simulation(20*10**6, demography)
print(get_migrating_tracts(ts))

Thanks for your help,

Leonardo

[jeromekelleher]

Hi @leovallini, welcome 👋

I haven't run your code or thought too hard about it, but your approach basically seems right to me. I think it'd be helpful if you printed out the migration events that you get here, and a simulated tree (ts.draw_text() or ts.draw_svg()) (maybe for 2 samples each, from the populations). That would help us see what's going on.

[leovallini]

@jeromekelleher below is the list of migrations with dest = 2 (the id of population "C") and tree drawn using ts.first().draw_svg(). If you wish to reproduce the results I used random_seed = 17
In the lines in bold the node of the migration is not one of the sample nodes of the admixed population. In this case since there are just 2 samples in the admixed population, that means that both of them have that segment (which migrated from pop "C"), but how do I know which individuals of the admixed population have it if I have, say, 50 individuals in that population?

Migration(left=8363514.0, right=8434505.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=15)
Migration(left=8445205.0, right=8509755.0, node=37, source=3, dest=2, time=1000.0, metadata=b'', id=16)
Migration(left=4815208.0, right=4835846.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=17)
Migration(left=5141491.0, right=5344206.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=18)
Migration(left=5378108.0, right=5539249.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=41)
Migration(left=6080849.0, right=6081860.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=42)
Migration(left=1385663.0, right=1483755.0, node=41, source=3, dest=2, time=1000.0, metadata=b'', id=43)
Migration(left=1483755.0, right=1605910.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=44)
Migration(left=1610059.0, right=1702641.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=45)
Migration(left=1713648.0, right=1725403.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=46)
Migration(left=17579549.0, right=17778726.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=47)
Migration(left=3814020.0, right=4007553.0, node=34, source=3, dest=2, time=1000.0, metadata=b'', id=48)
Migration(left=19713127.0, right=19720711.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=67)
Migration(left=19757407.0, right=19921131.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=68)
Migration(left=19921131.0, right=20000000.0, node=24, source=3, dest=2, time=1000.0, metadata=b'', id=69)
Migration(left=16254574.0, right=16262338.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=70)
Migration(left=16262338.0, right=16282084.0, node=98, source=3, dest=2, time=1000.0, metadata=b'', id=71)
Migration(left=16282084.0, right=16316389.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=72)
Migration(left=19529469.0, right=19570586.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=73)
Migration(left=19689890.0, right=19713127.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=74)
Migration(left=12329856.0, right=12399637.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=76)
Migration(left=13119144.0, right=13195006.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=77)
Migration(left=18889077.0, right=19080567.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=92)
Migration(left=11810581.0, right=11873464.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=100)
Migration(left=12921736.0, right=12926725.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=101)
Migration(left=14701772.0, right=14708526.0, node=12, source=3, dest=2, time=1000.0, metadata=b'', id=103)
Migration(left=14091862.0, right=14125673.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=105)
Migration(left=17133244.0, right=17143883.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=106)
Migration(left=17761844.0, right=18020450.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=107)
Migration(left=18233900.0, right=18315340.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=158)
Migration(left=8899431.0, right=9112335.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=185)
Migration(left=11322702.0, right=11367352.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=186)
Migration(left=2289921.0, right=2348596.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=191)
Migration(left=564611.0, right=596080.0, node=14, source=3, dest=2, time=1000.0, metadata=b'', id=200)
Migration(left=596080.0, right=608590.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=201)
Migration(left=7814384.0, right=7840873.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=204)
Migration(left=11783070.0, right=11810581.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=205)
Migration(left=14943026.0, right=15046290.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=206)
Migration(left=10194417.0, right=10260918.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=210)
Migration(left=10889584.0, right=11022002.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=211)
Migration(left=11348043.0, right=11372859.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=212)
Migration(left=11041748.0, right=11093413.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=213)
Migration(left=9801560.0, right=9819503.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=220)
Migration(left=5946054.0, right=5971591.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=223)
Migration(left=12776818.0, right=12800116.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=225)
Migration(left=3020923.0, right=3175110.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=227)
Migration(left=11954317.0, right=12053450.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=232)
Migration(left=11773103.0, right=11783070.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=233)
Migration(left=12800116.0, right=12863376.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=234)
Migration(left=13522236.0, right=13557691.0, node=8, source=3, dest=2, time=1000.0, metadata=b'', id=235)
Migration(left=19720711.0, right=19757407.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=237)
Migration(left=1483755.0, right=1610059.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=242)
Migration(left=11013878.0, right=11041748.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=243)
Migration(left=8710007.0, right=8890993.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=244)
Migration(left=10642269.0, right=10730673.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=245)
Migration(left=10260918.0, right=10333053.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=260)
Migration(left=8760307.0, right=8881935.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=261)
Migration(left=9112335.0, right=9664325.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=262)
Migration(left=27362.0, right=63588.0, node=16, source=3, dest=2, time=1000.0, metadata=b'', id=268)
Migration(left=63588.0, right=165713.0, node=42, source=3, dest=2, time=1000.0, metadata=b'', id=269)
Migration(left=9668584.0, right=9746802.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=286)
Migration(left=4473293.0, right=5028684.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=291)
Migration(left=7703622.0, right=7814384.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=292)
Migration(left=8326291.0, right=8344291.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=293)
Migration(left=596080.0, right=605225.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=295)
Migration(left=5154096.0, right=5161183.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=296)
Migration(left=7080440.0, right=7221140.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=297)
Migration(left=7551085.0, right=7627491.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=306)
Migration(left=8238797.0, right=8252468.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=329)
Migration(left=9819503.0, right=9871796.0, node=0, source=3, dest=2, time=1000.0, metadata=b'', id=332)
Migration(left=15672210.0, right=15761697.0, node=1, source=3, dest=2, time=1000.0, metadata=b'', id=346)

[jeromekelleher]

I think the issue here is that your looking at the internal nodes of the tree - once the samples within a population have coalesced, they'll produce internal nodes, so you can't just match against the same IDs.

This would be easier to see if there wasn't so much recombination. What do you get if you turn off recombination entirely?

[leovallini]

with no recombination I guess individuals within the admixed population will migrate (backward in time) entirely (= all of their genome) either to population B (with prob=0.75) or to population C (with prob=0.25).
However if I turn off recombination I will not be able to reach my final aim: identifying within each single genome in the admixed population which segments come from which source population (you can think about it like a local ancestry problem) ... while writing that I had an idea: do you think looping through the internal nodes involved in migration with node.children() until I reach a sample node could do it?

[jeromekelleher]

I understand that you won't be able to answer the question you're ultimately interested in @leovallini, I was just suggesting this as a path towards understanding what's going on and developing a solution. I usually find it best to start with a single tree to understand what's going on, and then generalising.

What do you get if you print out the tree and migrations in this case?

[leovallini]

ok sorry, I hadn't understood that.
these are all the migrations I had and the tree when there is no recombination

Migration(left=0.0, right=20000000.0, node=7, source=3, dest=1, time=1000.0, metadata=b'', id=0)
Migration(left=0.0, right=20000000.0, node=2, source=1, dest=0, time=2000.0, metadata=b'', id=1)
Migration(left=0.0, right=20000000.0, node=8, source=1, dest=0, time=2000.0, metadata=b'', id=2)
Migration(left=0.0, right=20000000.0, node=6, source=2, dest=0, time=2000.0, metadata=b'', id=3)

[jeromekelleher]

Great! Now, it's probably a good idea to add in a few more samples so we get more coalescence. Also, printing out the node table will be handy here print(ts.tables.nodes) so that we can see which population each node comes from. You can also print out the migration table directly too.

[leovallini]

I am now simulating the following: samples = {"D": 50, "B": 10, "C": 10}

Nodes Table

id	flags	population	individual	time	metadata
0	1	3	0	0	b''
1	1	3	1	0	b''
2	1	3	2	0	b''
3	1	3	3	0	b''
4	1	3	4	0	b''
5	1	3	5	0	b''
6	1	3	6	0	b''
7	1	3	7	0	b''
8	1	3	8	0	b''
9	1	3	9	0	b''
10	1	3	10	0	b''
11	1	3	11	0	b''
12	1	3	12	0	b''
13	1	3	13	0	b''
14	1	3	14	0	b''
15	1	3	15	0	b''
16	1	3	16	0	b''
17	1	3	17	0	b''
18	1	3	18	0	b''
19	1	3	19	0	b''
20	1	3	20	0	b''
21	1	3	21	0	b''
22	1	3	22	0	b''
23	1	3	23	0	b''
24	1	3	24	0	b''
25	1	3	25	0	b''
26	1	3	26	0	b''
27	1	3	27	0	b''
28	1	3	28	0	b''
29	1	3	29	0	b''
30	1	3	30	0	b''
31	1	3	31	0	b''
32	1	3	32	0	b''
33	1	3	33	0	b''
34	1	3	34	0	b''
35	1	3	35	0	b''
36	1	3	36	0	b''
37	1	3	37	0	b''
38	1	3	38	0	b''
39	1	3	39	0	b''
40	1	3	40	0	b''
41	1	3	41	0	b''
42	1	3	42	0	b''
43	1	3	43	0	b''
44	1	3	44	0	b''
45	1	3	45	0	b''
46	1	3	46	0	b''
47	1	3	47	0	b''
48	1	3	48	0	b''
49	1	3	49	0	b''
50	1	1	50	0	b''
51	1	1	51	0	b''
52	1	1	52	0	b''
53	1	1	53	0	b''
54	1	1	54	0	b''
55	1	1	55	0	b''
56	1	1	56	0	b''
57	1	1	57	0	b''
58	1	1	58	0	b''
59	1	1	59	0	b''
60	1	2	60	0	b''
61	1	2	61	0	b''
62	1	2	62	0	b''
63	1	2	63	0	b''
64	1	2	64	0	b''
65	1	2	65	0	b''
66	1	2	66	0	b''
67	1	2	67	0	b''
68	1	2	68	0	b''
69	1	2	69	0	b''
70	0	3	-1	3.08682357	b''
71	0	2	-1	7.77452729	b''
72	0	3	-1	12.08910407	b''
73	0	2	-1	15.40963226	b''
74	0	3	-1	15.68411335	b''
75	0	3	-1	16.80680327	b''
76	0	3	-1	20.68316667	b''
77	0	3	-1	23.08003833	b''
78	0	3	-1	27.43741557	b''
79	0	3	-1	29.61446552	b''
80	0	3	-1	30.56444345	b''
81	0	3	-1	39.22583788	b''
82	0	3	-1	50.63617297	b''
83	0	3	-1	51.07805694	b''
84	0	3	-1	59.11054323	b''
85	0	1	-1	75.93873248	b''
86	0	3	-1	76.46725834	b''
87	0	3	-1	76.87483386	b''
88	0	2	-1	77.93954107	b''
89	0	1	-1	94.22001666	b''
90	0	3	-1	94.46867641	b''
91	0	3	-1	121.582389	b''
92	0	3	-1	122.9058323	b''
93	0	1	-1	124.05796612	b''
94	0	3	-1	155.17173516	b''
95	0	2	-1	167.60069479	b''
96	0	1	-1	169.86079225	b''
97	0	3	-1	181.68952289	b''
98	0	3	-1	184.7893674	b''
99	0	3	-1	188.67445377	b''
100	0	3	-1	214.54257414	b''
101	0	3	-1	220.2750083	b''
102	0	3	-1	224.60625265	b''
103	0	2	-1	237.8615159	b''
104	0	3	-1	254.92958745	b''
105	0	3	-1	271.98198087	b''
106	0	3	-1	288.65319199	b''
107	0	3	-1	309.72254043	b''
108	0	1	-1	326.43176502	b''
109	0	3	-1	329.60171632	b''
110	0	3	-1	343.6237467	b''
111	0	3	-1	347.77427155	b''
112	0	3	-1	351.92905379	b''
113	0	3	-1	358.00241587	b''
114	0	3	-1	387.2583558	b''
115	0	3	-1	428.9184023	b''
116	0	1	-1	470.70219913	b''
117	0	3	-1	473.5144356	b''
118	0	3	-1	523.3961659	b''
119	0	3	-1	663.97928483	b''
120	0	3	-1	682.9769813	b''
121	0	3	-1	693.57215997	b''
122	0	3	-1	737.25942625	b''
123	0	1	-1	1009.32928625	b''
124	0	1	-1	1027.18982548	b''
125	0	2	-1	1063.53696363	b''
126	0	1	-1	1099.08326624	b''
127	0	2	-1	1213.77429796	b''
128	0	1	-1	1233.0357156	b''
129	0	1	-1	1241.25280208	b''
130	0	2	-1	1249.39319314	b''
131	0	1	-1	1380.19751381	b''
132	0	1	-1	1528.34199683	b''
133	0	1	-1	1886.43283242	b''
134	0	0	-1	2154.76438793	b''
135	0	0	-1	2689.44074082	b''
136	0	0	-1	5636.62352342	b''
137	0	0	-1	6359.58913338	b''
138	0	0	-1	10126.72565843	b''

Migration Table

id	left	right	node	source	dest	time	metadata
0	0	20000000	122	3	1	1000	b''
1	0	20000000	97	3	1	1000	b''
2	0	20000000	105	3	2	1000	b''
3	0	20000000	114	3	1	1000	b''
4	0	20000000	113	3	1	1000	b''
5	0	20000000	121	3	1	1000	b''
6	0	20000000	119	3	1	1000	b''
7	0	20000000	120	3	1	1000	b''
8	0	20000000	133	1	0	2000	b''
9	0	20000000	114	1	0	2000	b''
10	0	20000000	132	1	0	2000	b''
11	0	20000000	125	2	0	2000	b''
12	0	20000000	63	2	0	2000	b''
13	0	20000000	130	2	0	2000	b''

[jeromekelleher]

Great! So, you're interested in the migrations with source=3 I think. Can you see where these nodes are on the tree, and how they relate to the samples you're interested in? Note the "population" column in the node table - this tells you the population that a node was "born" in.

[leovallini]

Thanks, I think I finally got it.
All I need to do is identifying all the sample nodes downstream of internal nodes which are present in the migration table.
Then when adding recombination I have to look at the tree[s] describing the genomic interval that migrated and for each save the coordinates and the sample nodes downstream of it.

[jeromekelleher]

Excellent! I think that because you have no migration in the model before the admixture you can just look directly at the migrating segments, but it's good to look at the trees as well just to confirm that this the right thing to do.

[leovallini]

This is one more run where I restored recombination, I used samples = {"D": 5, "B": 2, "C": 2}, and random_seed = 17,.

This is one of the migrations I get, where the node is an internal node and not one of the sample nodes, then information about the node and the 5 trees describing the corresponding interval:

Migration(left=14252017.0, right=14270154.0, node=242, source=3, dest=2, time=1000.0, metadata=b'', id=872)

Node(id=242, flags=0, population=3, individual=-1, time=898.57880882, metadata= b'')

If I got this correctly the migrating chunk with coordinates (left=14252017.0, right=14270154.0) existed 898.57 generations ago but was later broken down by recombination (at least in some individuals) into two smaller chunks: one with coordinates (left=14252017.0, right=14261228.0 - the first two trees) and one with coordinates (left=14261228.0, right=14270154.0 - the last three trees). The first of the smaller segments is now present in three of my sampled individuals from pop D (with id 0, 2 and 3) while the second one can be found in four individuals: 0, 2, 4 and 3.

So I should mark the segment(left=14252017.0, right=14270154.0) as beeing "inherited from pop C" in individuals 0, 2 and 3.
while individual 4 just has the second part of the interval(left=14261228.0, right=14270154.0).