Tskit: Support for reference sequences

[jeromekelleher]

For some applications it would be useful to know the reference sequence that a tree sequence coordinate space refers to. For example, with real data, we should (at a minimum) record the reference build (e.g., GRCh38) and the contig ID (e.g., chr22) associated with a tree sequence. Ideally, we would also like to be able to do things like:

for site in ts.sites():
     print(site.position, site.ref_allele)

As well as situations in which a well-known canonical reference is available, we may also have a one-off reference sequence that we wish to record, e.g., in simulations.

To support this, I suggest adding a reference section to the file store, with some fields. Roughly, these might looks like:

reference/build         -- e.g. GRCh38
reference/contig       -- e.g. chr22
reference/id              -- md5 hash of the sequence, c.f. refget
reference/sequence -- Actual sequence information

The references section would be itself optional (keeping backward compatability), and some fields would probably be optional within this section (for example, reference/sequence would definitely be optional --- no point in storing GRCh38 chr1 over and over again). We can imagine having some mechanisms for automatically retrieving sequences using refget, but this isn't at all necessary for a basic implementation.

At a high-level, we should try to follow any upstream standards as closely as possible, e.g. GA4GH refget and any others that are relevant.

Any thoughts @tskit-dev/all, @bhaller?

[hyanwong]

This seems very sensible, but I wonder if the word "reference" is a bit too context-dependent. Perhaps "genome_reference" or something similar would be clearer?

[bhaller]

In SLiM I'm calling it the "ancestral sequence". I'm fine with "reference" or "genome_reference" or any such variant. The rest of the proposal seems good to me. It should be clarified, when this is documented, exactly what format reference/sequence is expected to be in – FASTA? Is a FASTA header allowed? Are newline characters allowed (r, n)? Nucleotides only, presumably, not FASTA amino acid sequences? But FASTA codes for ambiguous nucleotides of various sorts? But I think SLiM will just dump an unbroken stream of ACGT into it, so that is all that I need to be supported.

[jeromekelleher]

Thanks @bhaller. To clarify, it'll be a flat array of chars, such that reference[pos] is the reference allele at position pos.

[petrelharp]

I think this is a great idea, and lean towards calling it just reference (although reference_sequence would also be fine). I don't think it should be called ancestral_sequence, because I'd argue these shouldn't be the same thing, necessarily. Here's some reasons:

• If we say they're supposed to be the same, then the ancestral_state column of the Site table is supposed to match the positions in the reference sequence, which is annoying and redundant.
• The reference sequence is usually not the ancestral sequence, in practice. (eg GRCh38), but is still useful.
• The reference sequence could still be used by e.g. SLiM to store, essentially, the ancestral sequence.

[hyanwong]

FWIW I would (marginally) vote for reference_sequence. reference sounds like it might contain metadata about bibliographic references relevant to the provenance of the data. But maybe I'm being too naive.

[molpopgen]

@petrelharp is right here. I actually get papers across my desk where someone has calculated Fay and Wu's H w.r.to the reference allelic state :(.

I do like @hyanwong's idea to call it reference_sequence. At first, I thought this thread was about adding bibtex support to the docs.

[petrelharp]

@bhaller is going to actually start storing a reference sequence on output from SLiM nucleotide models, so we should at least get it nailed down what we will call it. The proposal is, I think,

• reference_sequence/build -- e.g. GRCh38
• reference_sequence/contig -- e.g. chr22
• reference_sequence/id -- md5 hash of the sequence, c.f. refget
• reference_sequence/sequence -- Actual sequence information

All of these would be optional. How about it?

[jeromekelleher]

OK, sounds good. I haven't researched what the right names for build, contig etc should be, so perhaps we could just keep it to reference_sequence/sequence if we're going to ship files that include it?

Also reference_sequence/sequence seems repetitive and unintuitive. Maybe reference_sequence/data would be better?

[petrelharp]

Besides deciding what to name these entries in the kastore, we should record this somehow, possibly by providing API functions that will write and read it. We think the main question is: is the reference sequence part of a table collection object? It would seem natural to make it a part of the table collection, because then you would get it when you read the .trees file. But, in the SLiM application, the reference sequence is expected to be large, and making it part of the table collection object means that there will be more than one copy of it sitting around, because:

• we have to copy the table collection when we write it out; and also
• SLiM has an internal version of the reference sequence, so when reading in the .trees file with table_collection_load( ), if reference_sequence was part of the resulting object, it would need to be copied over to where SLiM wants it

Since a table_collection_t has a kastore_t, we could define a function tsk_table_collection_read_reference_sequence( ) that would go into the store and pull out the reference sequence (if indeed such a store is associated). I'm not sure how to deal with writing out, though?

[petrelharp]

We could at least write a function that would re-open a .trees file and add a reference sequence to it.

[jeromekelleher]

I don't know how we're going to manage it from the C perspective @petrelharp, and I'm afraid I just don't have time to work on it right now. I don't want to rush in changes to the C API without thinking through the consequences very carefully (particularly not with something this fundamental).

I think the best thing to do from a SLiM perspective is to reopen the kastore in append mode after tskit has written it out, and write the reference sequence/data key in. Should be ~10 lines of code.

[molpopgen]

Does it make sense to add just the reference state at each variable site vs
the entire reference?
__
Kevin Thornton
Associate Professor
Ecology and Evolutionary Biology
UC Irvine
www.molpopgen.org
github.com/molpopgen
github.com/ThorntonLab

On Tue, Mar 19, 2019 at 1:36 PM Peter Ralph [e-mail]
wrote:

We could at least write a function that would re-open a .trees file and
add a reference sequence to it.

—
You are receiving this because you are on a team that was mentioned.
Reply to this email directly, view it on GitHub
https://github.com/tskit-dev/tskit/issues/146#issuecomment-474572322,
or mute the thread
https://github.com/notifications/unsubscribe-auth/AGHnH_HtcVCCz56Nr_k_VW8NrK_Uscx8ks5vYUnIgaJpZM4bs1uU
.

[jeromekelleher]

Does it make sense to add just the reference state at each variable site vs the entire reference?

That was my original proposal, but Ben has reasons for wanting the full reference/ancestral sequence.

[bhaller]

I think the best thing to do from a SLiM perspective is to reopen the kastore in append mode after tskit has written it out, and write the reference sequence/data key in. Should be ~10 lines of code.

Yep, that's what we've got for now. I think Peter just wants to clean up the code; but that can wait indefinitely, what we have now is fine for release.

[petrelharp]

Sounds good. We're doing that; and can continue this discussion when (if) someone else needs it.

[petrelharp]

I'm redirecting discussion from #326 over here, where I think it belongs. There, @hyanwong wrote:

I think we need to distinguish between (a) tree sequences that using the (msprime default) infinite sites mutation model and (b) tree sequences that have sites at integer positions (and hence probably use finite mutation models ).

Instead of guessing based on whether positions are all integers we should instead do something cleaner. Here are two non-mutually-exclusive options:

1. have a reference_sequence property of tree sequences, which might be NULL, and then have things like haplotypes and fasta output optionally use this to output invariant sites also.
2. make a subclass of the tree sequence, say FiniteSitesTreeSequence (ick, though) that only has mutations (and also recombinations?) at integer positions

This issue is about (1); but it seems like we might want to do something like (2) to be able to do (1). Otherwise whenever we do something with the reference sequence we'll have to check if positions are all integers, and if not, throw an error (or, something?). At the python level everything would be cleanest, if the relevant columns were actually integer... but this would probably be a giant pain otherwise.

[hyanwong]

Instead of guessing based on whether positions are all integers we should instead do something cleaner.

Yes, that would be much better.

1. have a reference_sequence property of tree sequences, which might be NULL, and then have things like haplotypes and fasta output optionally use this to output invariant sites also.

Nice idea.

1. make a subclass of the tree sequence, say FiniteSitesTreeSequence (ick, though) that only has mutations (and also recombinations?) at integer positions

It's not strictly necessary to have recombinations at integer positions, of course, but I suppose it is the option most likely to be needed by users.

This issue is about (1); but it seems like we might want to do something like (2) to be able to do (1). Otherwise whenever we do something with the reference sequence we'll have to check if positions are all integers, and if not, throw an error (or, something?). At the python level everything would be cleanest, if the relevant columns were actually integer... but this would probably be a giant pain otherwise.

Are you suggesting that the position in a FiniteSitesTreeSequence is stored as a C double, but converted to an integer in the python accessor functions?

[hyanwong]

I think there is a difference between a finite sites tree sequence and a tree sequence with a reference genome, and the two issues are getting confounded here. It is perfectly possible to have a TS with sites at finite positions, but without any reference sequence available. Indeed, I'm working on hacked versions of such tree sequences at the moment.

I suggest we might want to open a separate issue about finite sites?