Selective censoring of a tree sequence

[LukeAndersonTrocme]

The ancestry documentation provides an example of how to remove the pedigree information from an output tree sequence as follows:

tables = ped_ts.dump_tables()
tables.individuals.clear()
tables.nodes.individual = np.full_like(tables.nodes.individual, tskit.NULL)
censored_ts = tables.tree_sequence()
censored_ts

I would like to remove the individual_names for each node, but keep other metadata (such as geographic coordinates, and marriage dates). I've tried censoring the tree sequence following the example, but it seems to delete all meta data.

Does it make more sense to selectively delete all instances of individual_name in the metadata?

In [1]: import msprime
   ...:import tskit
   ...:import numpy as np

In [2]: ts = tskit.load("test_with_geo_22_sim.ts")

In [3]: print("Metadata for individual 999:", ts.individual(999).metadata)
Metadata for individual 999: {'geo_coord': [-71.5799, 46.6983], 'individual_name': '1234567'}

In [4]: tables = ts.dump_tables()
   ...:tables.individuals.clear()
   ...:tables.nodes.individual = np.full_like(tables.nodes.individual, tskit.NULL)
   ...:censored_ts = tables.tree_sequence()

In [5]: print("Metadata for individual 999:", censored_ts.individual(999).metadata)
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-14-fba9f8a6cd4c> in <module>
----> 1 print("Metadata for individual 999:", censored_ts.individual(999).metadata)

/lustre06/project/6068353/python_environments/pedsim/lib/python3.9/site-packages/tskit/trees.py in individual(self, id_)
   5094             metadata,
   5095             nodes,
-> 5096         ) = self._ll_tree_sequence.get_individual(id_)
   5097         return Individual(
   5098             id=id_,

IndexError: record index out of bounds

CC: @sgravel @jeromekelleher

[edit: the good news is that I figured out how to add the geographic coordinates to each individual!]

[hyanwong]

This is a tskit problem isn't it, rather than an msprime one?

There's an example of editing one field of the metadata here:

https://tskit.dev/tutorials/tables_and_editing.html#minor-edits

[petrelharp]

Yep - that documentation says how to do this sort of thing. Let us know if you run into problems! (And - I might transfer this to a tskit discussion.)

[LukeAndersonTrocme]

Just reading through the docs now, this is exactly what I was looking for. Thanks for the quick reply @hyanwong !

I'll post a brief follow up once I've implemented the changes in case anyone else finds this thread in the future.

(Apologies for placing this issue in the wrong channel)

[petrelharp]

No problem! There's no wrong place for questions, it's easy for us to re-file them.

[hyanwong]

No problem. Best to ask than not to ask at all.

[LukeAndersonTrocme]

Thanks for all the help! The docs show a clear example for modifying metadata. I've done essentially the same thing, but just for the sake of pasting an answer to the initial problem:

    tables = ts.dump_tables()
    new_metadata = [{"individual_name": "censored",
                     "geo_coord":i.metadata["geo_coord"]} for i in tables.individuals]
    validated_metadata = [
        tables.individuals.metadata_schema.validate_and_encode_row(row) for row in new_metadata
    ]
    tables.individuals.packset_metadata(validated_metadata)
    
    censored_ts = tables.tree_sequence()

    print("Metadata for individual 999:", censored_ts.individual(999).metadata)
Metadata for individual 999: {'geo_coord': [-71.5799, 46.6983], 'individual_name': 'censored'}

Two questions:

1. Is there a cleaner way of just removing one field of metadata but keeping the rest?
2. How can I remove the parents for each individual? (Since we've already removed the identifying IDs, I don't think this is absolutely necessary, but just to be overly cautious it might be worth doing.)

[LukeAndersonTrocme]

hmm that last one doesn't seem to be working. Am I doing something wrong?

def del_individual_name(md):
    del md["individual_name"]
    return md
    
def censor_pedigree(ts):
    tables = ts.dump_tables()
    new_metadata = [del_individual_name(i.metadata) for i in tables.individuals]
    validated_metadata = [
        tables.individuals.metadata_schema.validate_and_encode_row(row) for row in new_metadata
    ]
    tables.individuals.packset_metadata(validated_metadata)
    # remove parents
    tables.individuals.packset_parents=([] * tables.individuals.num_rows)
    censored_ts = tables.tree_sequence()
    return(censored_ts)
    
In [36]: print("Metadata for individual 500:", censored_ts.individual(500).metadata, "parents for individual 500:", censored_ts.individual(500).parents)
Metadata for individual 500: {'geo_coord': [-71.0233, 46.9703]} parents for individual 500: [ 86 151]

In [37]: print("Metadata for individual 500:", ts.individual(500).metadata, "parents for individual 500:", ts.individual(500).parents)
Metadata for individual 500: {'geo_coord': [-71.0233, 46.9703], 'individual_name': '1234567'} parents for individual 500: [ 86 151]    

[hyanwong]

Oh, parents is an array of arrays, right? so maybe

packset_parents=([[]] * tables.individuals.num_rows)

or

packset_parents=([np.array([], dtype=int)] * tables.individuals.num_rows)

[LukeAndersonTrocme]

Yeah so I think packset_parents is doing what we expect it to do, but I think the issue might lie in second line where we try to replace the values.

packset_parents=([np.array([], dtype=int)] * tables.individuals.num_rows)
tables.individuals.packset_parents=packset_parents

In [23]: packset_parents[1:5]
Out[23]: 
[array([], dtype=int64),
 array([], dtype=int64),
 array([], dtype=int64),
 array([], dtype=int64)]

[hyanwong]

Doh, my fault. This was a typo - there shouldn't be an =: it's a function

tables.individuals.packset_parents([[]] * tables.individuals.num_rows)

[LukeAndersonTrocme]

This works! Thanks for taking the time to walk me through this