tools for multiple chromosome

[petrelharp]

Once SLiM moves to the 'one ts per chromosome' model, we'll need some new tools. The main consideration here is that the nodes/indivs/pops tables need to agree between all the chromosomes' tree sequences.

For starting a SLiM simulation with msprime, I see two options:

1. simulate the whole genome in msprime then chop it up into a bunch of separate tree sequences
2. simulate each chromosome independently then union together the nodes tables, effectively saying that the sets of ancestors are disjoint between chromosomes

[hyanwong]

I am thinking about this in the light of some questions asked about multiple chromosomes on the SLiM mailing list. I think the main thing we might need is a way to record the lengths of the chromosomes somewhere (e.g. in metadata). We could then have a function to read in a "trees archive" as a single (concatenated) tree sequence, storing the sequence lengths somewhere, then another function to cut up the concatenated tree sequence back into separate ones of the sort loaded from one of the .trees files in the SLiM archive. Here's something I worked up:

import tskit
import json
from itertools import accumulate
import os

import numpy as np


def _extract_this_chromosome(ts):
    return json.dumps(ts.metadata['SLiM']['this_chromosome']), ts

def ts_from_archive(archive_name):
    chroms = dict(
        _extract_this_chromosome(tskit.load(os.path.join(archive_name, filename)))
        for filename in os.listdir(archive_name)
        if filename.endswith(".trees")
    )
    # Check all the chromosomes have identical metadata before merging
    global_metadata = None
    nodes_table = None
    individuals_table = None
    populations_table = None
    for ts in chroms.values():
        md = ts.metadata
        if "SLiM" not in md:
            raise ValueError()
        if "chromosomes" not in md["SLiM"] or 'this_chromosome' not in md["SLiM"]: 
            raise ValueError()
            
        del md['SLiM']['this_chromosome']
        if global_metadata is None:
            global_metadata = md
            nodes_table = ts.tables.nodes
            individuals_table = ts.tables.individuals
            populations_table = ts.tables.populations
        else:
            if md != global_metadata:
                raise ValueError
            if ((ts.tables.nodes != nodes_table) or
                (ts.tables.individuals != individuals_table) or
                (ts.tables.populations != populations_table)
            ):
                raise ValueError
    # Put in the same order as in the metadata
    chroms = [chroms[json.dumps(k)] for k in global_metadata['SLiM']["chromosomes"]]
    # Add the sequence lengths to the metadata
    seq_lens = [ts.sequence_length for ts in chroms]
    for d, l in zip(global_metadata['SLiM']["chromosomes"], seq_lens):
        d['sequence_length'] = l
    tables = chroms[0].dump_tables()
    # metadata for the final ts should not have a 'this_chromosome'
    schema = ts.metadata_schema.schema
    schema['properties']["SLiM"]['required'].remove('this_chromosome')
    schema['properties']["SLiM"]['properties']['chromosomes']['items']['required'].append("sequence_length")
    tables.metadata_schema = tskit.MetadataSchema(schema)
    tables.metadata = global_metadata
    ts = tables.tree_sequence()
    node_map = np.arange(ts.num_nodes)
    return ts.concatenate(
        *chroms[1:],
        add_populations=False,
        node_mappings=[node_map]*(len(chroms)-1),
        record_provenance=False)


def split_ts(ts, **kwargs):
    # take a concatenated tree sequence from SLiM and split it into multiple ones
    # keyed by the metadata['SLiM']['this_chromosome']['symbol'] value
    # **kwargs are passed to simplify()
    schema = ts.metadata_schema.schema
    schema['properties']["SLiM"]['required'].append('this_chromosome')
    schema['properties']["SLiM"]['properties']['chromosomes']['items']['required'].remove("sequence_length")
    md = ts.metadata
    lengths = [chrom.pop('sequence_length') for chrom in md["SLiM"]["chromosomes"]]
    breaks = list(accumulate([0] + lengths))
    tree_seqs = {}
    for lft, rgt, chr_md in zip(breaks[:-1], breaks[1:], md["SLiM"]["chromosomes"]):
        tables = ts.keep_intervals([(lft, rgt)], simplify=False, record_provenance=False).dump_tables()
        tables.trim(record_provenance=False) # NB - may not work as expected if chromosome has missing regions
        tables.simplify(**kwargs, record_provenance=False)
        tables.metadata_schema = tskit.MetadataSchema(schema)
        md['SLiM']['this_chromosome'] = chr_md
        tables.metadata = md
        tree_seqs[chr_md['symbol']] = tables.tree_sequence()
    return tree_seqs

So you can use it like:

# Yan_trees_archive from https://github.com/tskit-dev/tskit/issues/176#issuecomment-2959763856
merged_ts = ts_from_archive("/Users/yan/Downloads/Yan_trees_archive")
tree_seqs = split_ts(merged_ts)

print(tree_seqs["2"])

For the moment I have stored the chromosome lengths in e.g. metadata["SLiM"]["chromosomes"][0]["sequence_length"], but It might be good to have a non-SLiM-specific place to store chromosome lengths if we want to generalise this to other tools.

[bhaller]

Once SLiM moves to the 'one ts per chromosome' model, we'll need some new tools. The main consideration here is that the nodes/indivs/pops tables need to agree between all the chromosomes' tree sequences.

For starting a SLiM simulation with msprime, I see two options:

1. simulate the whole genome in msprime then chop it up into a bunch of separate tree sequences
2. simulate each chromosome independently then union together the nodes tables, effectively saying that the sets of ancestors are disjoint between chromosomes

For (1) does it work in msprime to specify a rate map with r=0.5 points to create independently segregating pseudo-chromosomes? If so, maybe it would be nice to have a utility function of some kind that sets up that structure, simulates it, and then chops it up for you at those positions?

Annotation of sets of tree sequences, comprising a trees archive, with the appropriate metadata like the chromosomes and this_chromosome key should also be a part of this work, presumably.

[bhaller]

I am thinking about this in the light of some questions asked about multiple chromosomes on the SLiM mailing list. I think the main thing we might need is a way to record the lengths of the chromosomes somewhere (e.g. in metadata). We could then have a function to read in a "trees archive" as a single (concatenated) tree sequence, storing the sequence lengths somewhere, then another function to cut up the concatenated tree sequence back into separate ones of the sort loaded from one of the .trees files in the SLiM archive. Here's something I worked up:

...

Something like this certainly seems useful. I can't really comment on the specifics of this implementation since I'm not really much of a user of the Python end of things. The thread on slim-discuss that I recall was about recapitation in particular, so it'd be good to figure out how recapitation fits into this picture (in terms of how the joined tree sequence gets recapitated easily, what recombination rate map it should use, etc.). An example vignette or two would be good.

For the moment I have stored the chromosome lengths in e.g. metadata["SLiM"]["chromosomes"][0]["sequence_length"], but It might be good to have a non-SLiM-specific place to store chromosome lengths if we want to generalise this to other tools.

I guess there are two questions here. One is where to put the sequence lengths used by your code when it joins/splits. That could go somewhere else in top-level metadata, certainly. Since this is a facility specifically for the benefit of your join/split mechanism, maybe it should go inside its own top-level metadata key? e.g., metadata["MultiChromJoinSplit"]["sequence_lengths"]?

The other question is whether the multi-chromosome metadata more generally would be better placed outside of SLiM's top-level metadata – not just the lengths, but the index, id, symbol, and type metadata that SLiM stores, too. It's all in SLiM metadata at present because at the time this work was happening in SLiM there seemed to be little time/interest in having any kind of support for multiple chromosomes in tskit, but the mood seems to have shifted somewhat! :-> If the metadata scheme being used by SLiM seems good, we could simply move that metadata to a different location in future (and add the sequence length, if desired). (Perhaps saving it to both places for some period of time, for backward compatibility, but deprecating the location inside the SLiM metadata; or perhaps just making a clean break, if we think that few external dependencies on this have formed in the short time since SLiM 5.0 came out.) I would like to see the concept of a "trees archive" become the official tskit stance on multi-chromosome support, and then we could have the tskit doc say "Here's how a trees archive should be configured, here's the metadata used, etc.". That would be best for inter-operability of tools in a multi-chromosome world. None of this implies that tskit itself would add any built-in support for trees archives or operating on multiple chromosomes; but declaring a standard for it would still be helpful as a point of reference for everybody. SLiM and pyslim can't really make such a declaration, it seems to me, since we're both in the SLiM world; we're a spoke in the wheel, the hub of which is tskit is. Standards should be declared at the level of the hub, even if the hub doesn't implement the standard in question. Or so it seems to me. :->

[hyanwong]

Yes, I agree, having a standardised way of representing multiple chromosomes is really a tskit thing. Having said which, I'm in favour of trialling the routines in SLiM/pyslim, and when the dust settles, implement compatible functionality in tskit.

For efficiency reasons, it would be nice (although not strictly necessary) to allow multiple tree sequences on disk to share the same node, population, etc tables. This would be a nice guarantee that they could be merged. I'm not sure quite how this would work, though, and is something for the distant future, I suspect.

[bhaller]

Yes, I agree, having a standardised way of representing multiple chromosomes is really a tskit thing. Having said which, I'm in favour of trialling the routines in SLiM/pyslim, and when the dust settles, implement compatible functionality in tskit.

OK, sure. But it might be good to make this clear to other folks, perhaps, to avoid everybody striking of in their own different directions? We want to avoid having the way that multiple chromosomes are handled becoming Balkanized.

For efficiency reasons, it would be nice (although not strictly necessary) to allow multiple tree sequences on disk to share the same node, population, etc tables. This would be a nice guarantee that they could be merged. I'm not sure quite how this would work, though, and is something for the distant future, I suspect.

Yes. For now, there is a lot of duplicated data, but gzip seems to compress it out to some extent?

Regarding what you wrote on slim-discuss about the correct rate for msprime for free recombination between chromosomes being log(2) rather than 0.5, interesting. The same would be true in SLiM also, except that we do some math behind the scenes to compensate; you can read about that in the SLiM manual in section 14.11, starting perhaps three pages into that section with "This is as good a place as any for a digression that I think this manual ought to contain somewhere". Just FYI. :->

[gregorgorjanc]

We will need multiple chromosomes soon and would indeed like to avoid Balkanisation! I believe Nemo devs are thinking about this too. Should we form a working group to this effect?

[hyanwong]

My ideal would be to store the information in the top tree sequence metadata for now, with a view to eventually make it a top-level attribute, like time_units. I'm not clear in my head how this works with things like sex chromosomes though, so I think experimenting via metadata is probably worth doing first.

I believe Nemo devs are thinking about this too

Is Nemo using the tskit format then?

[bhaller]

I believe Nemo devs are thinking about this too

Is Nemo using the tskit format then?

They are planning on doing so; there was a meeting between me, Peter, Gregor's postdoc Jinyang, and Fred Guillaume (the Nemo author) just recently.

[hyanwong]

I just put up a message on the "general" slack channel. I think if would be helpful to have somewhere to open multiple issues: that could be this GitHub repo, or we could create a new (temporary) one so that we can test-drive code there?

I think a good start would be to make a list of functions we might need to convert between the "multiple tree sequences in an archive" format, and the "single concatenated genome" approach. There's a lot of useful discussion in tskit-dev/tskit#176. I also recall that there was a way in the C layer for multiple table collections to share the same table (e.g. the nodes table): it would be nice to be able to create, save, and load immutable tree sequences with those guarantees, but I don't know how it would be done. Meanwhile I think that simply comparing tables to make sure they are equal would be fine, performance-wise (which is what I do in the code above)