htsget: spec should indicate (dictate?) whether .vcf.gz is preferred (mandatory?) over plaintext .vcf?

[brainstorm]

As a result of this discussion, the htsget-rs team at @umccr gave some thought to it and decided to favor indexed (and compressed) variant formats with htsget (.vcf.gz{.tbi}) for at least a couple of reasons:

1. Sets a bad precedent w.r.t the use of compression in production. It should be a given that assets are stored compressed at rest: it just makes economic sense at scale.
2. Supporting plaintext VCF unnecessarily complicates resolution logic on htsget server(s).

I'm sure there are more pros, cons, opinions, corner cases, workarounds and preferences about this topic?

EDIT: Tangentially related issue on how to handle ID resolvers and the requested format's file extensions (or lack thereof): umccr/htsget-rs#127

/cc @jrobinso @mmalenic @ohofmann @andrewpatto

[jrobinso]

I see, but this possibility is not mentioned in the spec, it probably should be. It was unexpected to request VCF format and receive a binary response.

Dealing with bgzip as a client can be complicated. Is it to be expected that the response(s) once concatenated contain a series of valid and complete gzip blocks?

WRT expectations for someone who is not deeply involved, I would expect the server to do any decompression needed and return the response in the format requested. I would even expect it to trim features (variants, alignments) that fall outside the query range requested. This is not possible with bgzip, of course, you can only return entire blocks. That would be some "value add" for using an htsget server, as it is this is going to be at least as complex and just accessing the file directly using the index.

The storage format on the server is a logically separate issue, and of no real interest to clients anyway. I would certainly expect VCF to be stored compressed.

[jrobinso]

Also, just a hint, when returning a request for a header don't assume the header is entirely contained in the first bgzip block. It will sometimes extend into the second block, or even third in extreme cases. Speaking from experience here.

[andrewpatto]

I always thought the returning VCF bgzip blocks was the "variants" mirror of returning compressed BAM blocks for "reads".

So the expectation is that if possible the htsget server could choose to select blocks by byte range out of the large object in its object store (I mean totally it could also choose to rewrite content and serve it up altered).

[andrewpatto]

So yes - I think it is definitely expected that the response(s) once concatenated would contain a series of valid and complete blocks (as in - I don't think the client has to be any cleverer than the BAM case)

[jrobinso]

@andrewpatto That's fine.. I suggest this be added to the spec though, there is no mention of it that I can find. It does complicate things for clients, one of the potential value adds of using an htsget server vs just accessing the file directly with an index would be that it deals with all this.

[andrewpatto]

Yes sorry - not arguing that the spec shouldn't be clearer on this! Just reflecting on how I had thought about it..

[jrobinso]

Yes it should be clearer, or at least mentioned as a possibility.

I can implement support for this even without the spec, but is IGV the only client? For that matter does UMCCR have the only server?

[jrobinso]

@brainstorm Please open issues for IGV and igv.js, if you care, to support bgzipped VCF from htsget servers. Supporting this in IGV is likely to be more complicated than igv.js, it probably will need to involve the htsjdk.

[brainstorm]

Please open issues for IGV and igv.js

Thanks Jim, will do!

I would expect the server to do any decompression needed and return the response in the format requested.

Well, BAM, CRAM, BCF are all compressed, so why do you have that particular expectation for VCF? Also, doing on-the-fly decompression is quite detrimental to bandwidth, latency and potentially overall system perf in general?

I would even expect it to trim features (variants, alignments) that fall outside the query range requested. This is not possible with bgzip, of course, you can only return entire blocks. That would be some "value add" for using an htsget server(...)

Agree, that'd be part of #838 most probably.

[jrobinso]

@brainstorm I have that expectation from the spec. Also, all previous servers that I am aware of were implemented that way. I do think your server is out-of-spec, but this might be a semantic distinction. Also, AFAIK you have the only htsget server, and are the only users of it, so in a practical sense whatever you do is the spec.

[jkbonfield]

Just a note, tbi isn't ideal as (I think) it limits chromosomes to 512MB. So requiring it may be detrimental in some circumstances. Bcftools supports CSI too which overcomes this limitation.

It depends whether you view htsget to be human only. IMO that would be a mistake.

Besides, what difference does the backend implementation make to the specification? The spec is the API access, not the format at rest.

[jmarshall]

Annoyingly, there's no explicit discussion in the relevant PRs (#233 and #301) and meeting minutes from 2017/18/19 of whether we were envisaging plain text VCF or bgzipped VCF being transmitted over the wire. My personal expectation (/ vague recollection perhaps) is that we meant plain text VCF, as that's the format named in the spec, and perhaps as that would be an MVP initially. During development of the original PR, we added BCF as an option but sadly there appears to have been no explicit discussion of compressed-VCF. It may be worth listening to some of the meeting recordings of the time, and would certainly be worth looking at original implementation behaviour, in particular of the DNAnexus implementation (if that supports/-ed the variants endpoint).

The htsget protocol was designed to be agnostic of how servers would store their data, in particular allowing both database-backed and file-backed information stores. As htsget fundamentally encourages genomic-region-based queries, the data storage would need to facilitate region queries. For file-backed servers (probably the gross majority!), that means that the underlying BAM/CRAM/VCF/BCF files will be compressed and indexed. But of course what appears on the wire in the protocol is a separate concern, and any indexing is an implementation detail that does not appear in the protocol.

Dealing with bgzip as a client can be complicated. Is it to be expected that the response(s) once concatenated contain a series of valid and complete gzip blocks?

Fundamentally an htsget client retrieves data blocks from the URLs listed in the ticket and concatenates them “to obtain the complete response in the requested format”. That is, it's extremely fundamental to the protocol that the concatenated response will be a valid file of the requested format. It's not stated explicitly, but this obviously implies that (for a bgzipped format) the concatenated response will consist of valid and complete BGZF blocks. This is most obvious in the context of requesting BAM format.

I would expect the server to do any decompression needed and return the response in the format requested.

An htsget server is indeed expected to do any transcoding implied by the request. For reads, it may need to convert between BAM and CRAM to satisfy the request's desired format. Or if it only has BAM present in its underlying data store, it might return an error for a request for CRAM. And vice versa. (I always assumed that generous servers might store their data in both BAM and CRAM to be able to satisfy both sorts of requests without transcoding each time.)

Similarly for variants, a file-based server will have indexed compressed-VCF and/or BCF as its backing file store and may need to transcode between them (and also to plain text VCF) to convert the internal data to the client's requested format. Unfortunately this doesn't answer the question of what exactly the client is asking for (plain text or bgzipped) when they request format=VCF!

BAM, CRAM, BCF are all compressed, so why do you have that particular expectation for VCF? Also, doing on-the-fly decompression is quite detrimental to bandwidth, latency and potentially overall system perf in general?

BAM, CRAM, and BCF are the protocol values the client can request with format=BAM|CRAM|BCF, and they are also the names of a bunch of formats that are unambiguously compressed binary formats. On the other hand, VCF — like SAM — is the name of an uncompressed plain text format — albeit one that is usually actually stored in compressed form in a .vcf.gz file. So I would share @jrobinso's expectation that format=VCF would be requesting the uncompressed plain text.

This is of course inefficient all around, hence my suspicion that the variants endpoint remains an uncompleted minimal-viable-product. It is true that all this deserves clarification and probably the addition of a way for clients to request bgzipped-VCF.

I would even expect it to trim features (variants, alignments) that fall outside the query range requested. This is not possible with bgzip, of course, you can only return entire blocks.

Yes, that is exactly what the requested query range is for — to omit reads/variants well outside the requested query range. In the BAM file-backed context, the anticipated implementation is that the server will use a bai/csi index to identify the BGZF blocks that contain reads within the query range and return a ticket corresponding to those blocks. (Recall that (sufficiently short) reads are not broken between BGZF blocks in all but the most ancient BAMs.) The htsget guarantee is that the BAM content returned contains the reads in the region you asked for and perhaps some additional reads too. If a server wanted to go the extra mile, it could decode the blocks on the edges of the region and return new BGZF blocks containing only reads that actually touch the query region — but this would be quite annoying to implement and is not currently expected.

Similar applies for servers backed by indexed bgzipped-VCF or BCF files. The csi/tbi index is used to identify the desired BGZF blocks, and the contents of those are transmitted back to the client. For format=BCF requests, the BGZF blocks can be returned just as for BAM. For format=VCF if that is plain text, the BGZF blocks would be decompressed and there would be an easy opportunity to return only the exact matching VCF lines. (See also #580 which would add a format=SAM mode which would similarly enable the server to return only the exact matching SAM records. Not sure what the ticket would need to look like though…)

(There is also samtools/htslib#1369 which prevents HTSlib from splitting VCF lines across BGZF blocks, which is an interesting counterpart to the longstanding not splitting BAM records across BGZF blocks. However this facilitates picking BGZF blocks based on an index for a compressed VCF file in the server backend, and again does not provide us with an answer re the format on the wire.)

Just a note, tbi isn't ideal as (I think) it limits chromosomes to 512MB

Like BAI, the TBI index format is fixed to a 14,5 geometry so is indeed limited to 512 MB. But as noted, none of the terms “bai”, “csi”, “tbi”, or “index” appear in the htsget spec or are exposed in the protocol. So the exact form of any indexing is merely a server implementation detail and is immaterial to the protocol.