ts: How will we encode node and individual IDs?

[gregorgorjanc]

Can we have just about any individual IDs in tskit tree sequences? I am asking because in AlphaSimR we have IDs as 1:n by default, growing forwards in time. Note though the pop@id and pop@iid (the later being integer) ... (I have not used this much so I am not 100% up to speed here.)

I think we will get 1:nSampledIndividuals from msprime by default, which could just be used in AlphaSimR as input and then we grow IDs from that point onwards. There is a glitch though where this approach will just work when we have nSampledIndividuals from the current population, but if we say use pedigree option in msprime then we will have 1:nSampledIndividuals oriented backwards in time and we will continue with nSampledIndividuals:(nSampledIndividuals+nAdditionalIndividuals) in AlphaSimR. I think this will all work well due to pop@id and pop@iid, but we should be mindful of this ...

Not sure what does SLiM do regarding individual IDs. @LynxJinyangii can you please check and report here?

@gaynorr can you advise what would be the best way forward here?

[LynxJinyangii]

Hi Gregor, I think SLiM just keeps its own stable pedigree IDs internally (https://benhaller.com/slim/SLiM_Manual.pdf, P448) and most of the information is stored in metadata (https://tskit.dev/pyslim/docs/latest/tutorial.html#individual-metadata). I am not entirely sure what your underlying concern is, but I think, aside from using metadata, we could also store this directly in the individual table. There is no requirement that the lengths of the individual table and node table match; the linkage between nodes and individuals is established via the 'individual' column in the node table, which references the corresponding row in the individual table. Not sure if this is what you asking, but a small example (if I use https://github.com/LynxJinyangii/AlphaSimRTmp/blob/main/dev/testData/_AlphaSimR_extended_chr1.trees):

ts1 = tskit.load('/Users/jliang2/R_scripts/AlphaSimRTmp/dev/testData/_AlphaSimR_extended_chr1.trees')
ts1.tables.nodes[491]

NodeTableRow(flags=0, time=-1.0, population=-1, individual=-1, metadata={'alphaSimR': {'id': '3_1'}})

print(ts1.tables.individuals)

╔══╤═════╤════════╤═══════╤═════════════════╗
║id│flags│location│parents│metadata         ║
╠══╪═════╪════════╪═══════╪═════════════════╣
║0 │    0│        │   2, 3│ {'file_id': '0'}║
║1 │    0│        │   4, 5│ {'file_id': '1'}║
║2 │    0│        │   6, 7│ {'file_id': '2'}║
║3 │    0│        │   8, 9│ {'file_id': '3'}║
║4 │    0│        │   6, 7│ {'file_id': '4'}║
║5 │    0│        │ 10, 11│ {'file_id': '5'}║
║6 │    0│        │ -1, -1│ {'file_id': '6'}║
║7 │    0│        │ -1, -1│ {'file_id': '7'}║
║8 │    0│        │ -1, -1│ {'file_id': '8'}║
║9 │    0│        │ -1, -1│ {'file_id': '9'}║
║10│    0│        │ -1, -1│{'file_id': '10'}║
║11│    0│        │ -1, -1│{'file_id': '11'}║
╚══╧═════╧════════╧═══════╧═════════════════╝

tables = ts1.dump_tables()
# add a new row at 12th for node 491 (ind ID:3; Hap: 1, in alphaSimR)
tables.individuals.add_row(parents = [2,1])
# assign this row to node 491
tables.nodes[491] = tables.nodes[491].replace(individual=12)
edited_ts = tables.tree_sequence()

edited_ts.tables.nodes[491]

NodeTableRow(flags=0, time=-1.0, population=-1, individual=12, metadata={'alphaSimR': {'id': '3_1'}})

edited_ts.tables.individuals[12]

IndividualTableRow(flags=0, location=array([], dtype=float64), parents=array([2, 1], dtype=int32), metadata={})

[gregorgorjanc]

@LynxJinyangii I see I write the issue in haste and without being clear on my concern. The concern I had/have is that we have ids oriented 1:n forwards in time (n is the youngest), while I see many applications in Python having them oriented 1:n backwards in time (n is the oldest). I imagine that SLiM has ids oriented 1:n forwards in time (n is the youngest) too somehow, but then I don’t know how ids are added together into one coherent set of ids during recapitation. You seem to be suggesting that we can have just about anything for ids and that ordering is done from topological sort, not numeric sort?

Could you create a small tree sequence with few individuals in a couple of different cases?

1. msprime: say 5 individual samples from the current population (these will have ids 0, 1, 2, 3, and 4, I believe) from msprime

2. msprime: say 5 individual samples from the current population (these will have ids these will have ids 0, 1, 2, 3, and 4, I believe) but also add 5 more based on a pedigree that connects back in time the 5 current pop samples (say, these could be parents or grandparents or such, and so these will have ids 5, 6, 7, 8, and 9, I believe)

3. SLiM: create 5 individual samples from the current generation (these will have ids these will have ids 0, 1, 2, 3, and 4, I believe) but also add 5 more based on a pedigree that connects forward in time from the 5 current pop samples (say, these could be children, and so these will have ids 5, 6, 7, 8, and 9, I believe)

You see what is starting to puzzle me?

4. Combine 3) with 1 or two, how does the tree sequence look like?

In fact I would like to understand how tree sequence will look like for each of these cases, in particular how we will handle ids going then forwards in time in AlphaSimR from either of the four cases above.

Maybe everything will just work out fine, I just haven’t done enough of tree sequence work in both directions to know, hence this issue to flesh this out;)

[LynxJinyangii]

@gregorgorjanc I think now I see what you mean. Just to clarify, individual IDs and node IDs can be completely independent. And, individual IDs are not mandatory in the tree sequences; they merely help us record pedigrees (for msprime, it only record that for the last two generations by default). So here, I assume what are you more concerned about are the node IDs. I think 1) the default behavior of these softwares is not the main point (in the tskit document, they state repeatedly that "The ordering requirements for a valid tree sequence do not specify that rows in the node table have to be sorted in any particular order, e.g. by time"); 2) we can easily adjust this with https://tskit.dev/tskit/docs/stable/python-api.html#tskit.TreeSequence.simplify based on our purposes. Here, I would like to keep using _AlphaSimR_extended_chr1.trees as an example:

1. this is the first tree in the tree sequence after adding generations by AlphaSimR:

print(ts1.first().draw_text())
                                                                       298                                                         
                                  ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┻━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓                   
                                 53                                                                            ┃                   
              ┏━━━━━━━━━━━━━━━━━━━┻━━━━━━━━━━━━━━━━━━━━┓                                                       ┃                   
             21                                       12                                                      19                   
              ┃                        ┏━━━━━━━━━━━━━━━┻━━━━━━━━━━━━━━━┓                                       ┃                   
              3                        2                               0                                       1                   
 ┏━━━┳━━━┳━━━┳┻━━━━┳━━━━━━━┓     ┏━━━┳━┻━┳━━━┓   ┏━━━┳━━━┳━━━━━━━━━┳━━━┻━━━━━━━━━┳━━━━━━━━━━━┓       ┏━━━┳━━━┳━┻━━━━━━━━━┓         
491 493 505 497   499     501   495 503 507 509 496 506 510       498           500         502     494 504 508         492        
             ┃   ┏━┻━┓   ┏━┻━┓                               ┏━━━┳━┻━┳━━━┓   ┏━━━╋━━━┓   ┏━━━╋━━━┓               ┏━━━┳━━━╋━━━┳━━━┓ 
            511 512 520 523 529                             513 515 517 519 514 516 518 521 525 527             522 524 526 528 530

2. add a line in individual table to see how it change (in case you're also interested in pedigree):

tables = ts1.dump_tables()
# add a new row at 12th for node 491 (ind ID:3; Hap: 1, in alphaSimR); note: I also added the metadata here for sorting
tables.individuals.add_row(parents = [2,1], metadata={'file_id':12})
# assign this row to node 491
tables.nodes[491] = tables.nodes[491].replace(individual=12)
edited_ts = tables.tree_sequence()
print(edited_ts.tables.individuals)
╔══╤═════╤════════╤═══════╤═════════════════╗
║id│flags│location│parents│metadata         ║
╠══╪═════╪════════╪═══════╪═════════════════╣
║0 │    0│        │   2, 3│ {'file_id': '0'}║
║1 │    0│        │   4, 5│ {'file_id': '1'}║
║2 │    0│        │   6, 7│ {'file_id': '2'}║
║3 │    0│        │   8, 9│ {'file_id': '3'}║
║4 │    0│        │   6, 7│ {'file_id': '4'}║
║5 │    0│        │ 10, 11│ {'file_id': '5'}║
║6 │    0│        │ -1, -1│ {'file_id': '6'}║
║7 │    0│        │ -1, -1│ {'file_id': '7'}║
║8 │    0│        │ -1, -1│ {'file_id': '8'}║
║9 │    0│        │ -1, -1│ {'file_id': '9'}║
║10│    0│        │ -1, -1│{'file_id': '10'}║
║11│    0│        │ -1, -1│{'file_id': '11'}║
║12│    0│        │   2, 1│  {'file_id': 12}║
╚══╧═════╧════════╧═══════╧═════════════════╝

3. if I want to order the nodes from the founder haps in AlphaSimR (use default params here):

# AlphaSimR generated nodes have ids from 491 to 530
edited_ts_simp1 = edited_ts.simplify(samples=list(range(491, 531)))
# now we can see that some of the nodes were removed because they don't have offspring in [491:530]
len(edited_ts.tables.nodes)
531
len(edited_ts_simp1.tables.nodes)
512
# And the rows of the removed nodes also disappeared in the individual table (and, at the last row, the parent IDs changed from [2,1] to [-1,1]!)
print(edited_ts_simp1.tables.individuals)
╔══╤═════╤════════╤═══════╤════════════════╗
║id│flags│location│parents│metadata        ║
╠══╪═════╪════════╪═══════╪════════════════╣
║0 │    0│        │ -1, -1│{'file_id': '0'}║
║1 │    0│        │ -1, -1│{'file_id': '1'}║
║2 │    0│        │ -1, -1│{'file_id': '6'}║
║3 │    0│        │  -1, 1│ {'file_id': 12}║
╚══╧═════╧════════╧═══════╧════════════════╝
# compare to ts1 and edited_ts, node 491 -> node 0
print(edited_ts_simp1.first().draw_text())
                                                  319                                        
                        ┏━━━━━━━━━━━━━━━━━━━━━━━━━━┻━━━━━━━━━━━━━━━━━━━━━━━━━━┓              
                       74                                                     ┃              
         ┏━━━━━━━━━━━━━━┻━━━━━━━━━━━━━┓                                       ┃              
         ┃                           44                                       ┃              
         ┃                ┏━━━━━━━━━━━┻━━━━━━━━━━━┓                           ┃              
        43               42                      40                          41              
┏━┳━━┳━━┳┻━━┳━━━━━┓   ┏━━┳┻━┳━━┓ ┏━━┳━━┳━━━━━━┳━━━┻━━━━━━┳━━━━━━━━┓    ┏━━┳━━┳┻━━━━━━━┓      
0 2 14  6   8    10   4 12 16 18 5 15 19      7          9       11    3 13 17        1      
        ┃  ┏┻━┓  ┏┻━┓                     ┏━━┳┻━┳━━┓  ┏━━╋━━┓  ┏━━╋━━┓          ┏━━┳━━╋━━┳━━┓
       20 21 29 32 38                    22 24 26 28 23 25 27 30 34 36         31 33 35 37 39
# the information inside just updated automatically, but the individual ID changed after simplification
edited_ts.tables.nodes[491]
NodeTableRow(flags=0, time=-1.0, population=-1, individual=12, metadata={'alphaSimR': {'id': '3_1'}})
edited_ts_simp1.tables.nodes[0]
NodeTableRow(flags=1, time=-1.0, population=-1, individual=3, metadata={'alphaSimR': {'id': '3_1'}})

4. if I want to order the nodes from the founder haps in AlphaSimR (use "filter_individuals=False" here):

edited_ts_simp2 = edited_ts.simplify(samples=list(range(491, 531)), filter_individuals=False)
len(edited_ts.tables.nodes)
531
len(edited_ts_simp2.tables.nodes)
512
# now some nodes are still removed, but won't cause any change in the individual table (Note: the IDs and parent IDs are independent from that in the node table.)
print(edited_ts_simp2.tables.individuals)
╔══╤═════╤════════╤═══════╤═════════════════╗
║id│flags│location│parents│metadata         ║
╠══╪═════╪════════╪═══════╪═════════════════╣
║0 │    0│        │   2, 3│ {'file_id': '0'}║
║1 │    0│        │   4, 5│ {'file_id': '1'}║
║2 │    0│        │   6, 7│ {'file_id': '2'}║
║3 │    0│        │   8, 9│ {'file_id': '3'}║
║4 │    0│        │   6, 7│ {'file_id': '4'}║
║5 │    0│        │ 10, 11│ {'file_id': '5'}║
║6 │    0│        │ -1, -1│ {'file_id': '6'}║
║7 │    0│        │ -1, -1│ {'file_id': '7'}║
║8 │    0│        │ -1, -1│ {'file_id': '8'}║
║9 │    0│        │ -1, -1│ {'file_id': '9'}║
║10│    0│        │ -1, -1│{'file_id': '10'}║
║11│    0│        │ -1, -1│{'file_id': '11'}║
║12│    0│        │   2, 1│  {'file_id': 12}║
╚══╧═════╧════════╧═══════╧═════════════════╝
# the individual ID won't change after simplification
edited_ts.tables.nodes[491]
NodeTableRow(flags=0, time=-1.0, population=-1, individual=12, metadata={'alphaSimR': {'id': '3_1'}})
edited_ts_simp2.tables.nodes[0]
NodeTableRow(flags=1, time=-1.0, population=-1, individual=12, metadata={'alphaSimR': {'id': '3_1'}})

5. if I want to order the nodes from the founder haps in AlphaSimR (use "filter_nodes=False" here):

# nothing will happen,
edited_ts_simp3 = edited_ts.simplify(samples=list(range(491, 531)), filter_nodes=False)
len(edited_ts_simp3.tables.nodes)
531
# even the order of the nodes won't change
print(edited_ts_simp3.first().draw_text())
                                                                       298                                                         
                                  ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┻━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓                   
                                 53                                                                            ┃                   
              ┏━━━━━━━━━━━━━━━━━━━┻━━━━━━━━━━━━━━━━━━━━┓                                                       ┃                   
              ┃                                       12                                                       ┃                   
              ┃                        ┏━━━━━━━━━━━━━━━┻━━━━━━━━━━━━━━━┓                                       ┃                   
              3                        2                               0                                       1                   
 ┏━━━┳━━━┳━━━┳┻━━━━┳━━━━━━━┓     ┏━━━┳━┻━┳━━━┓   ┏━━━┳━━━┳━━━━━━━━━┳━━━┻━━━━━━━━━┳━━━━━━━━━━━┓       ┏━━━┳━━━┳━┻━━━━━━━━━┓         
491 493 505 497   499     501   495 503 507 509 496 506 510       498           500         502     494 504 508         492        
             ┃   ┏━┻━┓   ┏━┻━┓                               ┏━━━┳━┻━┳━━━┓   ┏━━━╋━━━┓   ┏━━━╋━━━┓               ┏━━━┳━━━╋━━━┳━━━┓ 
            511 512 520 523 529                             513 515 517 519 514 516 518 521 525 527             522 524 526 528 530

So, I think one way to define the node order could be use the 'samples' list in simplify() (if we want them to decrease monotonically backward in time, we can use the last generation from AlphaSimR as samples). I am not sure if we can reoder the nodes without any deletion of them (maybe not necessary at the end), but if what we care about is the order of individual IDs/ pedigree (not that of node IDs), we can order/ fix it independently from the nodes. Another related function might be https://tskit.dev/tskit/docs/stable/python-api.html#tskit.TreeSequence.union

[gregorgorjanc]

@gregorgorjanc I think now I see what you mean.

Excellent;)

Just to clarify, individual IDs and node IDs can be completely independent.

Yes.

And, individual IDs are not mandatory in the tree sequences; they merely help us record pedigrees (for msprime, it only record that for the last two generations by default).

Yup.

So here, I assume what are you more concerned about are the node IDs.

I was more concerned about individuals (I always think at individual level first, but appreciate that ARGs are all about nodes first). But you raise an important point here! We could get individual IDs from msprime/SLiM, but we will always get node IDs from msprime/SLiM. This means that we must have a way of handling individual IDs and node IDs in AlphaSimR such that we can work with msprime & SLiM without issues. I have now added nodes to the title of this issue!

For example I am thinking of such cases:

• initiate ARG in AlphaSimR, maybe recapitate it with msprime, possibly add few more future generations with SLiM
• generate ARG with msprime, hand it over to AlphaSimR to add a few more future generations
• etc.

Can you think of a proposal on how we should handle these IDs to make this work smoothly?