2
$\begingroup$

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.

Improve this question
$\endgroup$
5
  • $\begingroup$ samtools and tabix support remote ranged file queries, so why bother implementing it yourself? $\endgroup$
    – Devon Ryan
    Jan 28, 2020 at 7:38
  • 1
    $\begingroup$ It supports remote ranged queries only for ftp: http: I picked a poor use case example. I actually need to access the file over a protocol not supported by samtools -- S3. It is https, but has custom authentication headers. I might also want to do the same over ssh or some other protocol, or via custom authentication scheme, like ntlm or whatnot. I was trying to avoid making this an s3 or protocol question -- I just need the file offsets. $\endgroup$
    – init_js
    Jan 28, 2020 at 8:44
  • $\begingroup$ Here is the full documentation: htslib.org/doc/htslib-s3-plugin.html $\endgroup$
    – Devon Ryan
    Jan 28, 2020 at 8:49
  • $\begingroup$ ugh. they're going to support every protocol under the sun eventually? There's already very good s3 clients out there. Will it support datasync, parallel uploads, separation of each gzip block into its own s3 object multipart fragment too to take advantage of better replication? Strange decision imo. What if my bam file's split in a million files? $\endgroup$
    – init_js
    Jan 28, 2020 at 8:51
  • $\begingroup$ It's supported via a plugin, so presumably there will just be more of those in the future. $\endgroup$
    – Devon Ryan
    Jan 28, 2020 at 8:53

1 Answer 1

Reset to default
2
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ From the htslib side of things there'd be no reason to load an index without a corresponding data file, since there are already functions for performing the range queries. $\endgroup$
    – Devon Ryan
    Jan 28, 2020 at 11:12
  • 1
    $\begingroup$ For htslib, the APIs you would be interested in are those discussed in htslib #385. For various reasons, something similar will likely make it into HTSlib in the not too distant future. $\endgroup$
    – John Marshall
    Jan 28, 2020 at 12:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.