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.
samtools
andtabix
support remote ranged file queries, so why bother implementing it yourself? $\endgroup$