Finalise semantics of fasta output

[jeromekelleher]

Pending some issues being worked out about what do when we have floating point site positions and integer site positions (see #326, #350)

[hyanwong]

To address this issue, I do think it would be helpful to know what we plan to do about integer site positions. Do we have a consensus here? E.g. we could make a derived class, e.g. FiniteSitesTreeSequence as suggested by @petrelharp. If so we could store the site positions as true integers, or simply as floating point (and report them as integers). Alternatively, we could keep the same class, and have a flag in the metadata that tells us to interpret the existing floating point positions as integers. Other suggestions welcome too.

I'm hoping to get @dianalir to have a go at thinking about this soon, so comments very welcome, @jeromekelleher

[jeromekelleher]

Integer site positions aren't necessary, I think.

I think the semantics are: if there is no reference sequence present (which there never is, right now; see #146 for discussion) we output the haplotypes as computed by the haplotypes() method. If there is a reference sequence, we use this to fill in the gaps between the variant sites.

But for now, I think it's clear: we just output the haplotypes.

[hyanwong]

OK. For the record, I think that it is quite misleading to concatenate together non-adjacent sites as adjacent letters in a FASTA file. The FASTA convention is (AFAIK) for letters to represent adjacent bases, and for potential gaps (of any length) to be represented by a -. So my eventual preference would be to output a set of non-adjacent sites something like:

A-T-T-C-G-... etc.

I do think that would be the least misleading and most standards-conformant output. Of course, once we have integer sites (and/or reference sequences: I don't think they are mutually dependent), then in some circumstances we could omit the - character between sites.

Perhaps I'm being too pernicketty, though?

[hyanwong]

(I forgot that I wrote up my thoughts on this in #326 (comment))

[jeromekelleher]

With non-integer site positions then there's no real meaning to 'between' them. We can add the padding behaviour later as an option, but I think the only reasonable thing to do for now is just output the haplotypes. It's up to the user to put in intermediate sites if they want them.

[jeromekelleher]

Hmmm, reading the fasta "specifications" available I'm coming around to @hyanwong's point of view. A "-" character does indicate a "gap of indeterminate length" in these, so it would be a good choice to (as Yan puts it) indicate that this is weird fasta data. I've definitely seen alignment programs that use "-" to indicate a gap of length exactly one though, so I wonder how well this is supported/interpreted.

It might be a good exercise to see what various libraries do when they are fed with fastas containing these characters. For example, what would Pysam and BioPython make of a FASTA string "-A-T-"?

[hyanwong]

It is common to allow both . and - for these sort of gaps, and the A2M format uses these to distinguish - as character to represent a missing state (i.e. used as an "alignment column"), and a period (.) to "represent insertion positions between alignment columns". The A2M docs also state that "The ... periods ... are only for human readability, and may be omitted."

So my suggestion is as follows (happy to change anything though)

The haplotypes() function is given an optional extra parameter, e.g. pad_sites, which, if true, inserts a single dot (not a dash) between sequence letters if the site positions are non-integer (and possibly if they are integer but are missing some integer positions). However, this parameter defaults to False, so that our current behaviour is unchanged. Since haplotypes() doesn't claim to be FASTA format, this is fine.

The write_fasta() function by default calls haplotypes(pad_sites=True), so that unlike haplotypes(), the default FASTA output does have dots between letters if the site positions are non-integer. We probably want the ability to turn this off, however - perhaps by having the function definition something like write_fasta(omit_site_padding=False).

Either way, we still indicate samples that have missing data at a particular site with a -, in both haplotypes() and write_fasta(). There is an argument that users of haplotypes() might want to turn this off, if they only want a single sample haplotype, and don't care about aligning. This is a bit of an edge-case, however, and could be incorporated later if indeed there is any call for it.

[hyanwong]

p.s. I'm happy to implement my suggestion above (or whatever modification of it we agree upon), if I get the go ahead, @jeromekelleher. I can do this today, tomorrow, or Fri.

[jeromekelleher]

Right, well this is looking like rather a lot more work than just quickly popping out a point release. I think we'll just have to undocument the fasta functionality for now, since it's controversial.

[hyanwong]

Yes, happy to work on this is there's agreement, with the understanding that it's pushed to 0.2.4. I think it's better to wait a little and get the format right, rather than end up having to support something that's not quite what we want.