Use multi iterators in the synced reader

[daviesrob]

Adds bcf_itr_regions() and tbx_itr_regions() interfaces as BCF and tabix equivalents to sam_itr_regions().

Uses the new interfaces to add multi-iterator support to the BCF/VCF synced reader. This makes it considerably faster when index jumping over a large number of regions by reducing the amount of work done looking up index entries. Note that this only affects regions supplied via bcf_sr_set_regions() which works using iterators over the input files. The similar bcf_sr_set_targets() function works as a filter and does not use iterators.

The locations supplied to bcf_sr_set_regions() can either be an in-memory list, or read from a file (which can itself be indexed with tabix). This change should enable multi-iterators for both of these cases.

There are a few changes to behaviour as a result of this change, which unfortunately can't be completely hidden from callers due to the visibility of the synced reader structures in the public headers.

• The bcf_sr_regions_t::start and bcf_sr_regions_t::end fields don't track the location of the current output records as accurately as before because the iterators no longer stop at the end of each listed region.
• bcf_sr_regions_t::regs may be non-NULL when reading regions from a tabix indexed file.
• More storage is needed when reading regions from a tabix indexed file, as regions are now handled in chromosome-sized chunks instead of individually.

Additionally this fixes a bug where setting the BCF_SR_REGIONS_OVERLAP option to 2 ("true variant overlap") with bcf_sr_set_opt() could result in records being returned out of order. This could happen if the variant in the record occurred after the end of the first region it overlapped, but also hit a later region.

Fixes samtools/bcftools#2557

[jkbonfield]

How does this

The bcf_sr_regions_t::start and bcf_sr_regions_t::end fields don't track the location of the current output records as accurately as before because the iterators no longer stop at the end of each listed region.

Interact with the frameshifts plugin here:
https://github.com/samtools/bcftools/blob/develop/plugins/frameshifts.c#L125-L145

It explicitly checks start and end for each region. The comment is a bit vague about tracking as accurately as before, but it sounds like this should now fail. I'll explore more in testing.

Edit: it's a different interface and is apparently not interacting with the changes, although this is very unclear initially.

[jkbonfield]

Mostly minor issues. Typos and queries. Quite a few of my problems with it stem from the original code. I needed to figure out what was going on in order to understand the changes, and that's far from clear due to the opaqueness and lack of adequate commenting.

Comparing against bcftools/plugins/frameshifts.c, and our conversation earlier:

Initially it was unclear why frameshifts isn't bugged given the change in behaviour of bcf_sr_regions_t::start/end. However I believe it's the difference between bcf_sr_next_line vs bcf_sr_regions_overlap. Although it's unclear how bcf_sr_regions_next also differs here. I think it falls into the same camp as bcf_sr_regions_overlap, which is that it doesn't use a multi-region iterator.

It looks like bcf_sr_next_line operates on file(s), while bcf_sr_regions_next operates on region lists instead (agnostic of associated files). The documentation doesn't adequately explain why it's beneficial when processing a single file to add it to an array of incoming synced-readers and iterate that way instead of using the manual approach. It probably should add a comment however stating the benefits of doing multi-region work when using the file based interfaces. (Plus maybe a comment on bcf_sr_regions_overlap as to why it's not using a file, which is a bit confusing given the naming convention.)

Also "bcf_sr_regions_next" claims to return the next interval, while "bcf_sr_regions_overlap" checks a specific location, presumably (but unstated) updating bcf_sr_regions_t::start/end itself if the queried region does infact overlap. This needs documenting, as bcf_sr_regions_next currently does.

That said, if bcf_sr_regions_next isn't updating start/end correctly in some circumstances, then those also need explicitly mentioning in the documentation too. Or rather I suspect it's setting them correctly, but the regions may not be the same ones given in the input file due to region-merging?

[jkbonfield]

Typo: "Create"

[daviesrob]

Typo fixed.

[jkbonfield]

Any chance we can fix the function declaration to use something other than int i. Looking at the review I couldn't figure out what was happening, because my instinct was this is a loop variable yet we only have seq for that. I later realised i is passed in as an argument!

Something far more descriptive would be useful, especially given it has no documentation or comment above the function. It doesn't need to be long or overly verbose, just not something most commonly associated with trivial local looping.

[daviesrob]

It's now reader_idx and some documentation has been added to the header.

[jkbonfield]

This still leaves reglist->reg as a copy to the old pointer. Does it matter if we don't completely copy reglist_in and have a shared pointer?

I believe not, as it's just a cached copy of regions->seq_names[regions->iseq] anyway, but thought I'd double check.

[daviesrob]

No, it shouldn't matter. The pointer needs to be valid when calling tbx_itr_regions() / bcf_itr_regions() but after that it isn't used. Ownership of the string pointed to is also expected to be elsewhere.

[jkbonfield]

How is it possible for reg[idx].start > reg[idx].end? This feels like corrupted data, as each region should always have end >= start.

Is this a typo for something else?

Edit: no - it's documented in region1_t as:

    hts_pos_t start, end;   // records are marked for skipping have start>end

This probably needs a comment as the code simply looks wrong.

[daviesrob]

Comment added.

[pd3]

I checked: this is required as a consequence of regions_merge() which can produce redundant regions. It is cheaper to mark as redundant (by making them invalid) then memmove the structures around.

[jkbonfield]

I cannot find a function named merge_regions.

After trial and error I found regions_merge instead. This and the other mention need renaming. It's not relating to changes in this PR, but I was attempting to figure out why we checked start > end and discovered these comments in the process.

[daviesrob]

Changed to regions_merge(), along with the other instances mentioned.

[jkbonfield]

Fix to regions_merge().

[jkbonfield]

If we're modifying internal function declarations, then should we also make them C specification compliant by removing the leading underscores?

[daviesrob]

Done.

[jkbonfield]

I still need to do some integration testing. The above was simply a review by reading the code, rather than stress testing it. That's up next. I'll put on my evil-user hat! :-)

[jkbonfield]

If we have a bed file that contains no overlap, eg our VCF is chr20 only and our bed file is for chr1 only, then there's a difference, albeit in warning only.

Dev prints up nothing, while this PR prints up:

[W::hts_itr_regions] Region 'chr1' specifies an unknown reference name. Continue anyway

This is incorrect, as chr1 is infact in the header contig lines.

[jkbonfield]

I see vcfstats.c uses bcf_sr_next_line() , and hence speeds up with the multi-region iterator. However it also has a do_indel_stats function which accesses start and end directly.

I thought I managed to get a difference in output, but it was just me forgetting to grep out the header lines.

Nonetheless, it is running considerably faster in this PR, implying that the multi-region iterator is in use. So I wonder in what scenarios the start / end are incorrect and whether I'm just getting lucky on this data set.

It may be that the field it's accessing is the underlying bed iterator rather than the region iterator. It's all a bit convoluted!

[jkbonfield]

The following is a problematic BED file as start == end, and start is zero based so it should always be less than end.

chr1    3336336 3336338                                                                                  
chr1    3336345 3336345                                                                                  
chr1    3336345 3336353                                                                                  

Nonetheless, a zero-sized region should not invalidate the next line, and our handling of this case is a regression.

for i in dev pr;do echo === $i ===;./bcftools.$i view -R _2.bed /nfs/users/nfs_j/jkb/scratch/data/dragen-novaseq.pcr-free.20x/HG002.novaseq.pcr-free.20x.vcf.gz|egrep -v '^#'|cut -f 1-5;done
=== dev ===
chr1	3336344	.	ATG	A
=== pr ===

Although chr1 3336345 3336353 would appear to be chr1:3336346-3336353 in 1-based coordinations, as this is a 2bp deletion it spans that region and should be found.

If I remove either of the first two lines, then the 3rd line now reports the entry in both tools.

It's marginal though, so maybe this is an off-by-one error lurking in there somewhere.

[jkbonfield]

The above bug is fixed by this:

diff --git a/synced_bcf_reader.c b/synced_bcf_reader.c
index 6b3f193f..a05ef0f3 100644
--- a/synced_bcf_reader.c
+++ b/synced_bcf_reader.c
@@ -936,9 +936,12 @@ static int multi_check_overlap(region_t *reg, int reader_idx,
     if ( midx >= reg->nregs )
         return 0;  // pos is after all regions
 
-    // Skip records starting inside previous region or
+    // Skip records starting inside previous (non-skipped) region or
     // ending before the current one starts
-    hts_pos_t prev_end = midx > 0 ? reg->regs[midx - 1].end : -1;
+    int pidx = midx-1;
+    while (pidx >= 0 && reg->regs[pidx].start > reg->regs[pidx].end)
+        pidx--;
+    hts_pos_t prev_end = pidx >= 0 ? reg->regs[pidx].end : -1;
     if ( beg <= prev_end || end < reg->regs[midx].start )
         return 0;

The problem is prev_end was accessing a record which had been skipped due to it being zero length. The general skip strategy (start>end, which it will be for zero length data due to the nature of BED) fixes it.

[daviesrob]

The bug fix has been included and the "unknown reference" messages have been suppressed.

The accessing of bcf_sr_regions_t::start and bcf_sr_regions_t::end in bcftools stats and bcftools plugin frameshifts shouldn't be a problem. In both cases the bcf_sr_regions_t structure is created directly, and is being used via bcf_sr_regions_overlap() rather than to make an iterator. That usage is unaffected by this change.

[pd3]

I tested on a 24 GB VCF (vcf.gz) and 282,615 regions (tabix-indexed) using the command.

bcftools query -f '%CHROM\t%POS\t%REF\t%ALT' file.vcf.gz -R regions.txt.gz

The new version is consistently 4.5x faster.

The only difference between the outputs is that the original outputs duplicate regions.

Thank you for this!