checkDuckDBGRanges() started failing on the Bioconductor build machines
(26 test failures, all "Attributes: < Component 'row.names': Modes:
character, numeric >") after GenomicRanges 1.65.2 rewrote
as.data.frame.GenomicRanges(): it used to default row.names to
names(x) when not supplied, and now hardcodes row.names = NULL
internally (an explicitly-passed row.names argument is silently
discarded too), putting names(x) in a new names column instead.
checkDuckDBGRanges()'s own df <- as.data.frame(expected) reference
assumed the old convention, so its row names stopped matching
as.data.frame(object)'s, which still carries real names as row names.
Not a DuckDBGRanges bug: confirmed by installing the exact upstream
commit and reproducing, then verifying as.data.frame,DuckDBGRanges-method
itself was unaffected, since it goes through DuckDBDataFrame's own
conversion rather than GenomicRanges's. Fixed by explicitly restoring
rownames(df) <- names(expected) and dropping the new names column
before comparing, which works against both old and new GenomicRanges.
The same fix landed in BiocDuckDB's copy of this helper, and in
BiocDuckDB::writeParquet.GenomicRanges() itself, which had the same
reliance on the old convention and so silently wrote positional numbers
instead of real feature names as the __name__ column whenever the
source GenomicRanges had real names.nearest(), precede(), follow(), and distanceToNearest() returned
all-NA (or an empty Hits) for any DuckDBGRanges that does not carry an
explicit keycol, which is the default for an object built from a plain file.
Three separate defects were stacked here; the first two masked the third.
.add_keycol_indices() treated the keycols slot as a set of literal key
values and inner-joined the table against them. A row-number-keyed frame
does not store key values there: it stores the c(NA, -n) sentinel
returned by set_row_number(). The join therefore matched nothing and
silently emptied the result rather than erroring. The row position is now
derived directly instead. Note that the frame's own row_number column
cannot be used as the index as-is, since a subset keeps the parent's
numbering (g[3:4] carries 3, 4 for what are now positions 1 and
2); a row-number-keyed object always iterates in ascending row_number
order, so ranking over that column yields the required 1..n position..build_nearest_single_result() used the collected index directly as a
subscript. Now that the index is computed in SQL it arrives as
integer64, and subsetting with an integer64 vector silently returns
NA for every position instead of indexing. (The Hits-building path was
unaffected only because it already coerced with as.integer().)nearest() lacked the is.na() guard on the joined subject index that
distanceToNearest() already had. The setup left-joins on seqnames, so a
query range on a seqname with no subject produces a row whose subject
bounds are NULL; because DuckDB's greatest() skips nulls rather than
propagating them, greatest(NULL, NULL, 0) is 0, so that row scored a
perfect distance, won its own minimum-distance filter, and yielded a hit
to a NULL subject. The directional precede()/follow() paths were
never affected, as their validity comparisons against a NULL bound
evaluate to NULL and drop the row.Every pre-existing test in this family supplied an explicit keycol, which is
why none of this was covered. Reported in part as item R-G2 of the
six-architect cross-engine parity review; the other two defects were found
while reproducing it. All four generics are now checked against
GenomicRanges for both keying modes.
psetdiff() handled only one of the four ways y[i] can overlap x[i],
as its own header comment conceded ("assumes y is fully contained overlap
scenario or edge-aligned"). A pair that does not overlap at all came back
wider than x, and a pair where y covers x came back with a
negative width, which is not merely a wrong number but unmaterializable
(each range must have a non-negative width). A y sitting strictly
inside x cannot be expressed as a single range, so base refuses it;
this returned nonsense instead. The four shapes now match base:
(y_end + 1, x_end) for a left-edge overlap, (x_start, y_start - 1) for
a right-edge one, a zero-width range at (x_start, x_start - 1) when y
covers x, an unchanged x when the two do not overlap, and an error when
y falls strictly inside.
Three further defects in the same method:
seqnames and strand compatibility were ignored entirely, so ranges
were subtracted from one another across different chromosomes. Base
leaves x[i] untouched whenever the pair cannot be compared.ignore.strand was accepted and then never used. It is now honoured, and
validated..build_DuckDBGRanges()'s default
coordinate sort, which silently permuted it: result[i] was not
x[i] minus y[i]. As with pgap(), the result is now ordered by the
pairing index.Pairing also moves off a bare row_number() over DuckDB's undefined scan
order onto .add_keycol_indices(), matching the pgap() rewrite. That is
separately tracked as review item R-G4 for punion()/pintersect(), but
correct pairing is a precondition for this method being correct at all.
Reported as item R-G1 of the six-architect cross-engine parity review. The
pre-existing test covered only left-edge overlaps, which is the single
shape the old code got right; all four shapes are now checked against
GenomicRanges for both keying modes and both ignore.strand settings.
punion() and pintersect() shared the same defects psetdiff() had,
through their common .parallel_set_op() engine, and base treats the two
differently in a way that was not modelled at all:
pintersect() on a pair that does not overlap produced a negative
width, which is unmaterializable (each range must have a non-negative width). Base returns a zero-width range at (start(x), start(x) - 1),
which is now the clamp applied whenever the coordinate expressions would
invert.punion() silently spanned a gap between the two members of a pair.
Base errors unless fill.gap = TRUE.seqnames and strand compatibility were not checked, so both
operations combined ranges across different chromosomes. Base refuses
such a pair outright in punion(), and returns the zero-width range in
pintersect().ignore.strand was accepted and never used, as were pintersect()'s
strict.strand (which tightens compatibility so * no longer matches
+/-) and drop.nohit.ranges (which removes the zero-width results).
All three are now honoured and validated.result[i] was not x[i] combined
with y[i]. They are now ordered by the pairing index.distance(x, y) returned a number instead of NA for a pair on different
seqnames whenever either strand was *. The strand | chain was not
parenthesized inside the &, and dbplyr does not add the parentheses
itself, so under SQL's precedence the rendered predicate reduced to
"same seqname AND strands equal OR x is * OR y is *" and any * strand
made the pair valid on its own. The pre-existing test missed it because it
never combined a seqname mismatch with a * strand. This is the same
hazard already documented on pgap().
The remaining element-wise operations (punion(), pintersect(), and
distance()) now pair x[i] with y[i] through each object's recorded
key rather than a bare row_number() over DuckDB's undefined scan order,
completing the correction the pgap() rewrite began and psetdiff()
continued. Reported as item R-G4 of the six-architect cross-engine parity
review; the correctness defects above were found while fixing it. All
three are now checked against GenomicRanges with and without an explicit
keycol.
narrow() silently accepted requests base refuses, and in two of them
returned a range wider than the one it was given, which is the opposite
of what the function is for.
width fixes whichever side was left NA, so exactly one of
start/end must be NA. narrow(x, width = w) on its own had nothing
to anchor to and quietly anchored at the start, and supplying all three
of start/end/width quietly ignored end. Both now error, as in
base.narrow() may only shrink a range (base solves the request with
allow.nonnarrowing = FALSE), and may not invert one. Requests such as
narrow(x, start = 5, end = 200) on a 101-wide range previously returned
104-299, extending well past the original. These now error. The check is
per row, matching base, so it also catches a request that is legal for
some ranges in the object but not others; it costs one aggregate query.The review reported this as item R-G3, described as narrow(width = w)
leaving the coordinates unchanged while writing an inconsistent width.
That part does not reproduce: .modify_DuckDBGRanges_datacols() recomputes
end as start + width - 1 whenever only the width changes, so the object
stays self-consistent. The real divergence is the missing validation above.
Every solvable form of the call already agreed with GenomicRanges and
still does.
case_when is now declared in the @importFrom tags of trim(),
restrict(), reduce(), and setdiff(), which all use it. They had been
relying on the import declared by psetdiff(), which no longer needs it.coverage() for DuckDBGRanges, computed lazily via the same
delta-event + window-function cumsum() sweep-line pattern already used by
disjoin()/gaps(): only the compact per-seqname breakpoint table is
collected (bounded by the number of ranges, never by genome length), and a
real SimpleRleList is built from it, matching
GenomicRanges::coverage()'s value type (integer unless a non-integer
weight is used). shift and width (per-seqlevel, truncating or
zero-padding) are supported, including clipping a range that lands (partly
or entirely) before position 1 after shift, exactly as
GenomicRanges::coverage() does, rather than erroring. weight is
scalar-only for now (a per-range vector or mcols-column-name weight is not
supported, and errors clearly rather than risking a silent row-order
misalignment).coverage() on a DuckDBGRanges object crashed with a
C stack overflow: with no coverage() method of its own, it inherited
coverage,GenomicRanges-method, which calls
split(ranges(x), seqnames(x)), and ranges,DuckDBGRanges-method returns
a DuckDBDataFrame rather than real IRanges. The new method above takes
priority via S4 dispatch, so this no longer happens.pgap(x="DuckDBGRanges", y="DuckDBGRanges") is rewritten as a genuine, lazy
SQL implementation (no collect()). Previously the method built a SQL
pipeline for the non-overlapping-gap case, discarded it unused, and always
fell back to materializing both inputs and delegating to
IRanges::pgap(); every call, including through the DuckDBGRanges,GRanges
and GRanges,DuckDBGRanges methods, paid for a wasted SQL round trip and
never stayed lazy despite otherwise participating in the package's SQL
pushdown design. Also fixes a latent return-type inconsistency where a
zero-length call returned a DuckDBGRanges while every other call returned
a GRanges.
start = min(end(x), end(y)) + 1,
end = max(max(start(x), start(y)) - 1, min(end(x), end(y)))) matches
IRanges:::pgap,IntegerRanges,IntegerRanges, including its convention of
a zero-width range at the boundary for overlapping or adjacent pairs.x[i]/y[i] are paired by each object's own recorded key values via
.add_keycol_indices() (the same mechanism .setup_nearest_neighbor_join()
uses for nearest()/precede()/follow()), not a plain row_number();
a plain row_number() has no guaranteed correspondence to a table's
logical row order across separate query executions and can silently
mispair rows. Row-number-keyed (no explicit keycol) DuckDBGRanges
inputs now error clearly rather than risk a silently wrong pairing.x[i]/y[i] compatibility (matching seqnames, and unless
ignore.strand, compatible strand) is validated in SQL, erroring like
IRanges::pgap() on a mismatch; the new ignore.strand argument is
threaded through the DuckDBGRanges,GRanges/GRanges,DuckDBGRanges
methods as well.x/y pairing rather than
the lazy-builder helper's default coordinate sort, which would otherwise
silently permute the result.DuckDBGRanges, consistent with range() and
reduce(), rather than a materialized GRanges.pgap() tests for the DuckDBGRanges return type and added
coverage for the seqnames/strand compatibility errors (including
ignore.strand) and the row-number-keyed guard..-prefixed names should not be exported:
.keycols -> keycols, .has_row_number -> has_row_number,
.set_row_number -> set_row_number, and
.makePrettyCharacterMatrixForDisplay -> makePrettyCharacterMatrixForDisplay.
Updated the importFrom/exportMethods entries, the DuckDBGRanges and
DuckDBGRangesList method definitions, and all call sites. Requires the
renamed version of DuckDBDataFrame.sapply(..., simplify = FALSE) with setNames(lapply(...), ...) in
the mcols normalization (preserving the element names), and 1:h with
seq_len(h) in the display method, resolving the corresponding R CMD BiocCheck NOTEs.system2() instead of system() in inst/scripts/run_vignette_benchmarks.R.system.time() so the
runtimes are visible, and rewrote the summary section as prose.URL and BugReports fields to DESCRIPTION.?DuckDBGRanges).library(BiocStyle) call from the vignettes (the
BiocStyle::html_document output and ::-qualified helpers apply the style
without loading the package). Applied suite-wide for consistency with the
Bioconductor review of DuckDBDataFrame.\value sections (roxygen @return) to the DuckDBGRanges and
DuckDBGRangesList class and utility man pages, documenting the values
returned by the constructors, accessors, replacement, coercion, and
range-operation methods. Resolves the R CMD BiocCheck "missing \value"
WARNING..explainQuery()
returns DuckDB's plan for a lazy query, and a test uses it to assert the
overlap join stays an IEJoin / range join rather than degrading to a
NESTED_LOOP_JOIN over the full cross product (which OOMs on skewed inputs).
Interval overlap is a range join, not an ASOF (nearest-match) join. The
overlap join was factored into an internal .overlap_join_tbl() helper so its
plan can be inspected without collecting.precede()/follow() with select = "all" now return only the subject(s) at
the nearest distance (the ties), matching base GenomicRanges — previously
they returned every directional subject. The select = "first"/"last" paths
were already correct; the "all" branch was missing the min-distance filter.distance() now returns NA for a paired range on a different seqname or
(unless ignore.strand = TRUE) an incompatible strand (+ vs -; *
matches any strand), matching base GenomicRanges::distance. Previously it
ignored seqnames and strand and returned a bare coordinate gap. Added an
ignore.strand argument.narrow() now treats a negative start as a position counting back from
the range end (-1 is the last base), matching base IRanges::narrow
(solveUserSEW) and the existing negative-end behaviour. Previously a
negative start was applied as start + (start - 1), walking left of the
range start and producing an invalid interval. The start/end/width expressions
are rebuilt so every sign-and-width combination is correct; because the
datacols expressions flatten to SQL without parentheses, each is kept as a
column +/- scalar (never a subtraction of a compound). Oracle-tested against
base narrow across positive/negative start/end and their width
combinations.nearest(x) and distanceToNearest(x) (no subject) now exclude self-hits,
matching base GenomicRanges' drop.self = TRUE for the missing-subject form:
a range is never its own nearest neighbour, so its next-nearest is returned
instead of itself (distance 0). Previously both delegated to the two-argument
form against x itself, so every range matched itself. The explicit
nearest(x, x) / distanceToNearest(x, x) forms still keep self-hits, as in
base. (precede/follow were already correct — their strict inequalities
never admit a self-match.) Added a self-query test comparing to the base oracle.nearest/precede/follow/distanceToNearest now treat strand "*" as
compatible with any strand, matching base GenomicRanges: "+" pairs with
{+,*}, "-" with {-,*}, and "*" with all. The neighbour join used a
strict equi-join on strand, which silently dropped every "*" pair (a "*"
query against "+"/"-" subjects returned NA). The join is now on
seqnames only, followed by a strand-compatibility filter.precede() and follow() are now strand-directional. Base defines them in the
transcription direction — for a "-" strand query the roles of upstream and
downstream are reversed, and for a "*" query the direction is chosen per
subject strand. The previous implementation used a fixed genomic-coordinate
direction regardless of strand, so results were inverted on the "-" strand.
The convention is now selected per row (subj_start > end vs subj_end < start) via use_minus = strand == "-" | (strand == "*" & subj_strand == "-"),
ranking by the transcription-direction gap; ignore.strand = TRUE collapses to
the "+" convention as before. follow() keeps its base select = "last"
(largest index) tie-break, precede() its select = "first" (smallest index).
Added oracle tests against base GenomicRanges across +/-/* queries with
mixed-strand subjects, plus tie-break tests.resize(x, width, fix = "center") now matches base IRanges::resize():
new_start = start + (width(x) - width) %/% 2, anchored at the (strand-independent)
center with an exact integer width. The previous implementation used the midpoint
(start + end) / 2 with true division --- yielding fractional coordinates for
even-width ranges --- and center + width %/% 2 - 1, which returned width - 1
for odd target widths. Added fix = "center" cases to the resize tests.DuckDBGRangesList).GRanges on scATAC-seq (1M peaks) and variant (10M variants) scenarios,
rendered from precomputed results so the vignette builds quickly.LIST[] representation of grouped ranges,
for developers.inst/scripts/ with the offline benchmark generator
(run_vignette_benchmarks.R) and the vignette table helpers
(make_timings_table.R), following the HDF5Array performance-vignette
precompute pattern.