hts_itr_multi_t versus hts_itr_t

[jmarshall]

The documentation in htslib/hts.h of the multi-iterator functionality is minimal, so it is unclear how to use it and when exactly to use it rather than the existing hts_itr_t.

What is the reason for there to be a separate hts_itr_multi_t and a separate bunch of API functions for dealing with it? Is the functionality of hts_itr_t now a subset of that of hts_itr_multi_t, or does each flavour have things that the other can't do? I note that samtools/sam_view.c appears to have separate code for one region and multiple regions, which does not seem ideal.

Is the reason for the new separate hts_itr_multi_t struct simply a question of binary ABI compatibility and not wanting to add fields to the existing hts_itr_t?

If so, as HTSlib is currently in an ABI-breaking mode, it would be possible to unify the two forms of iterator by renaming hts_itr_multi_t to hts_itr_t, making hts_itr_multi_t a typedef alias of it, and adjusting the API function implementations accordingly, adding hts_itr_multi_next()/etc compatibility wrappers in htslib/hts.h as necessary.

Is this on the roadmap for the ABI-breaking HTSlib release?

[valeriuo]

Is the functionality of hts_itr_t now a subset of that of hts_itr_multi_t?

Pretty much so. hts_itr_t can be a bit faster when the target regions are sparsely spread out.

Is the reason for the new separate hts_itr_multi_t struct simply a question of binary ABI compatibility and not wanting to add fields to the existing hts_itr_t?

There is this aspect and the fact that, depending on the target regions, they give slightly different results. If the target regions overlap, hts_itr_t visits the overlapping section twice, while hts_itr_multi_t removes this overlap.

If so, as HTSlib is currently in an ABI-breaking mode, it would be possible to unify the two forms of iterator by renaming hts_itr_multi_t to hts_itr_t, making hts_itr_multi_t a typedef alias of it, and adjusting the API function implementations accordingly, adding hts_itr_multi_next()/etc compatibility wrappers in htslib/hts.h as necessary.

Is this on the roadmap for the ABI-breaking HTSlib release?

Not yet, but it is a fair candidate.

[jmarshall]

Thanks for the reply.

There is this aspect and the fact that, depending on the target regions, they give slightly different results. If the target regions overlap, hts_itr_t visits the overlapping section twice, while hts_itr_multi_t removes this overlap.

“Different results” here is presumably a question of how third-party code (or e.g. samtools view) uses the APIs, not anything inherent in their implementations. So samtools view foo.bam reg1 reg2 reg3, which I assume is what you're talking about here, could be implemented via hts_multi_itr_t using a loop

foreach region on the command line
    hts_multi_itr_t *itr = sam_itr_regions(idx, …, region[i], …);
    …run through itr…

or by making one hts_reglist_t containing all the regions and just having one hts_multi_itr_t iteration. It would be a bit silly to do the former, but this is how calling code would get reads in the overlapping regions multiple times or just once, according to the calling code's desired behaviour.

So if you renamed hts_multi_itr_t to hts_itr_t, the behaviour of existing loop-over-all-the-regions code would be unchanged, returning the overlapping sections multiple times.

Point is, users would probably prefer just one hts_itr_t that is more capable than to be tempted to write multiple separate code blocks to use hts_itr_t or hts_multi_itr_t (especially as the “sample code“ in samtools view appears to suggest doing just that). This ABI-breaking release is an opportunity to tidy these iterator types and functions up (while keeping compatibility wrappers/#defines for the _multi_ names).

hts_itr_t can be a bit faster when the target regions are sparsely spread out.

That would appear to be something that should be improved in the hts_multi_itr_t implementation 😄

[brentp]

I agree that the documentation and use cases are lacking for hts_itr_multi_t. For example, if I am doing 1000 queries spread relatively evenly over a 3Gigabase genome, would hts_itr_multi_t give any speed improvement?

[valeriuo]

I agree that the muti-region iterator hasn't been adequately documented nor advertised, as it was mainly created to deal with some performance limitation factors in samtools view, such as the inability to iterate over CRAM files and the unnecessary file seeks.
That being said, there are a few implementations you can use as a starting point, until the documentation of this API is brought to a decent standard.
The first is in "view" and the second is in "stats". There is also a test case in htslib.

If you put your queries in a BED file and run samtools view over your data, with and without the -M flag, you should get a rough idea about the speed improvement. If the queries are sorted and sufficiently spread out, so that they fall in different index entries, there will be no improvement.

[brentp]

ok. thanks for the clarification. So hts_multi_itr_t iterates over the full file and checks for overlap? vs hts_itr_t which uses index-skipping? If so, for me, that tidbit would be sufficient to add to the header file.

[daviesrob]

hts_multi_itr_t uses index skipping, but allows you to specify more than one location to visit. The result is more or less the same as using a filter on the file, so a read will only be output once even if it covers more than one region.

PR #718 now includes more documentation for the iterator functions. It will hopefully be merged early next week.

As we're in an ABI breaking mode at the moment, it's also possible that hts_multi_itr_t and hts_itr_t will get merged together into the same structure. If that works then it will be possible to use them pretty much interchangeably. One will become a typedef of the other; you'll get a single- or multi- region iterator depending on whether you call sam_itr_queryi, sam_itr_querys or sam_itr_regions; and sam_itr_next and sam_itr_multi_next will become the same function.

[jmarshall]

(:tada: in the hope that @daviesrob's last paragraph comes to pass — this is what this issue is actually about!)