How to obtain file offsets corresponding to range query in block-gzipped file

Question

How can we retrieve the file offsets of the zipped-chunks in a block-gzipped file that contain records for a region/window of interest (e.g. Chr01:100-Chr02:145)?

It seems like the CLI of samtools, or tabix should support that, but I can't find a way to just get the offsets in the datafile, without having the accompanying datafile at hand. I couldn't find a formal specification for those index files either.

Perhaps one can recommend a library that can parse those index files, in python, C, or golang? It's hard to accept that there isn't a reusable block-gzipped index parser somewhere, other than the one buried deep in samtools/htslib/gatk code.

edit

I'd initially given a small motivating example making use of file offsets to access a subportion of file on a remote server. Commenters were quick to point out that samtools supported a series of protocols for doing exactly that: ftp, http, and s3 (which requires https, which requires that samtools be compiled with libcurl support) Thanks! That example was detracting, so I removed it.

samtools view does fetch a subportion of a file (potentially from a remote location), but it does a few extra things, like decompressing the chunks retrieved. With this question, I'm interested in manipulating the gz chunks directly.

Answer 1

So if the purpose of obtaining file offsets is to access subportions of remote files via protocols that samtools supports, as mentioned in this comment, then you can make your life easier and use samtools/tabix directly on the cli, prefixing your datafile with the corresponding scheme.

One side effect of using the CLI however, is that the blocks retrieved will be decompressed (e.g. samtools view). You'd have to re-compress them if what you're trying to do is create a smaller "bam" file.

Keep in mind that if you retrieve all zipped chunks which cover your range of interest, then you might be extracting more data than you asked for. The first and last chunk may contain records outside your region, just because their neighboring records happened to be compressed in the same chunk. You'd have to be careful to filter those downstream.

Specification for the tabix format

If you absolutely need a custom solution, the tabix data structure format, authored by Heng Li, is on the samtools website:

https://samtools.github.io/hts-specs/tabix.pdf

A more in depth discussion of how the data structure works can be found here:

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3042176/

The tree-shaped binning index is a bit hard to understand from the text alone, but the tabix paper mentions that their approach is almost identical to the binning algorithm in this paper: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC186604/pdf/X6.pdf (page 1003). That paper has a nice illustration which makes things clearer.

Operation

Note: This is from my understanding of how the indices interoperate.

You don't need the data file to retrieve chunk offsets, you can do the work with just the tabix file. But it's complicated, and probably error prone if you're implementing this yourself.

Let's assume the query interval is Chr01:100-Chr02:2000.

Each sequence (Chr01, Chr02, etc) has its own pair of indices (binning index + linear index) in the tabix file. So you would break the query into two queries: Chr01:100-999999999999 and Chr02:0-2000, and concatenate the two results.

Within one sequence, to find all the blocks for open interval [X,Y), you would use the binning index to retrieve all the bins overlapping this interval (There is C code to do this in the data structure pdf). Then, you need to go through all the possible bins selected, and filter chunks. Not all chunks in a given bin will be overlapping with X, Y.

To filter chunks within a bin, you need the help of the linear index. There is one entry in the linear index per 16kbp of sequence positions. Each entry of the linear index contains the virtual offset of the leftmost record overlapping that 16kb window. I don't know what the entry contains if there is no record overlapping that 16kb window, presumably 0.

You would find the linear entry in which position X falls (X/16Kb), take that offset, call it O_x. You find the linear entry corresponding to position Y, and take that offset, call it O_y. When you go through the selected bins, filter out any chunk whose offset is < O_x, or whose offset >= O_y.

Assuming your datafile is sorted, the chunk range you're interested in is contained between the minimum offset and the maximum offset selected.

The virtual offsets are 64bit unsigned integers. The higher 48 bits of the virtual chunk offset corresponds to the file offset of the datafile (the other 16 you don't need. They are the offset of the record within that decompressed chunk).

Software

The software I have found typically mimics closely the hts_lib API. And hts_lib doesn't export a function which loads an index without also requiring a datafile. As far as I know.