Strange setdiff behavior - identified bug

[gevro]

Hi,
I'm seeing strange behavior in setdiff. GenomicRanges v1.58.0

Two GRangesLists, examine list element 16384 in each of them:

> x_list %>% class
[1] "CompressedGRangesList"
attr(,"package")
[1] "GenomicRanges"
> rem_list %>% class
[1] "CompressedGRangesList"
attr(,"package")
[1] "GenomicRanges"

> x_list[[16384]]
GRanges object with 1 range and 2 metadata columns:
      seqnames            ranges strand |
         <Rle>         <IRanges>  <Rle> |
  [1]     chr4 85143672-85143727      + |
                                          run_id        zm
                                        <factor> <integer>
  [1] Lib1 233832556
  -------
  seqinfo: 195 sequences (1 circular) from hg38 genome

> rem_list[[16384]]
GRanges object with 1 range and 3 metadata columns:
      seqnames            ranges strand |
         <Rle>         <IRanges>  <Rle> |
  [1]     chr4 85143672-85143727      + |
                                          run_id        zm       hit
                                        <factor> <integer> <logical>
  [1] Lib1 233832556      TRUE
  -------
  seqinfo: 195 sequences (1 circular) from hg38 genome

setdiff between the two GRangesLists, and examine list element 16384 in the result:

> out_list <- GenomicRanges::setdiff(x_list, rem_list, ignore.strand = TRUE)

> out_list[[16384]]
GRanges object with 1 range and 0 metadata columns:
      seqnames            ranges strand
         <Rle>         <IRanges>  <Rle>
  [1]     chr4 85143672-85143727      *
  -------
  seqinfo: 195 sequences (1 circular) from hg38 genome

The result in out_list[[16384]] should be empty.

Confirmation that seqnames were encoded identically:

> as.integer(seqnames(x_list[[16384]]))
[1] 4
> as.integer(seqnames(rem_list[[16384]]))
[1] 4

Direct setdiff between the individual elements works correctly:

> GenomicRanges::setdiff(x_list[[16384]], rem_list[[16384]], ignore.strand = TRUE)
GRanges object with 0 ranges and 0 metadata columns:
   seqnames    ranges strand
      <Rle> <IRanges>  <Rle>
  -------
  seqinfo: 195 sequences (1 circular) from hg38_SMaHT genome

There might be something fundamental I'm not understanding about setdiff. Or perhaps this is some weird bug.

[gevro]

Doing more investigation (with chatGPT's help), suggests this is indeed a bug in psetdiff.

> max(as.integer(seqnames(x_list@unlistData)))
[1] 23
> max(as.integer(seqnames(rem_list@unlistData)))
[1] 4

And then when deconstructGRLintoGR() creates temporary seqnames codes, it uses zero padding width that depends on only one of the input elements, so when the number of digits differ for the maximum encoded integer seqnames, then the temporary seqnames codes don't match.

That's my best guess at least.

Some more data that may be helpful:

> x_list %>% unlist %>% seqnames %>% as.integer %>% table
.
   1    2    3    4    5    7    8    9   10   11   12   13   14   15   16   17 
4808 2100  316  593 1328 5144 1464  914 1856 1042 6580  204  260  610  906 1088 
  18   20   23 
 124 1316  230 
> rem_list %>% unlist %>% seqnames %>% as.integer %>% table
.
4 
1

xgr  <- GenomicRanges:::deconstructGRLintoGR(x_list)
ygr  <- GenomicRanges:::deconstructGRLintoGR(rem_list)

> xgr[16384]
GRanges object with 1 range and 2 metadata columns:
      seqnames            ranges strand |
         <Rle>         <IRanges>  <Rle> |
  [1] 16384|04 85143672-85143727      + |
                                          run_id        zm
                                        <factor> <integer>
  [1] Lib1 233832556
  -------
  seqinfo: 30883 sequences from an unspecified genome


> ygr[1]
GRanges object with 1 range and 3 metadata columns:
      seqnames            ranges strand |
         <Rle>         <IRanges>  <Rle> |
  [1]  16384|4 85143672-85143727      + |
                                          run_id        zm       hit
                                        <factor> <integer> <logical>
  [1] Lib1 233832556      TRUE
  -------
  seqinfo: 1 sequence from an unspecified genome

You can see that the temporary seqnames is being created as 16384|04| from x_list and as 16384|4` from rem_list. The zero padding is inconsistent (that's the bug) and this is the causing setdiff to wrongly not subtract the ranges from each other.

Is there a workaround for this?

[gevro]

Reproducible simple example (with the help of chatGPT):

suppressPackageStartupMessages({
  library(GenomicRanges)
  library(GenomeInfoDb)
  library(S4Vectors)
})

## Shared Seqinfo with many seqlevels
si <- Seqinfo(
  seqnames   = paste0("chr", 1:23),
  seqlengths = rep(1e6, 23),
  genome     = "toy"
)

## x_list uses chr4 and chr23 (so max f2 = 23)
x_list <- GRangesList(
  GRanges("chr4",  IRanges(100, 120), strand = "+"),
  GRanges("chr23", IRanges(1,   1),   strand = "+")
)

## rem_list uses only chr4, and second element is empty (so max f2 = 4)
rem_list <- GRangesList(
  GRanges("chr4", IRanges(100, 120), strand = "+"),
  GRanges()
)

## IMPORTANT: make seqlevels identical BEFORE setting seqinfo
seqlevels(x_list)   <- seqlevels(si)
seqlevels(rem_list) <- seqlevels(si)

seqinfo(x_list)   <- si
seqinfo(rem_list) <- si

## Show the max integer code differs, even though seqlevels are identical
max(as.integer(seqnames(x_list@unlistData)))     # 23
max(as.integer(seqnames(rem_list@unlistData)))   # 4
identical(seqlevels(x_list), seqlevels(rem_list))# TRUE

## Expected: x_list[[1]] should be removed by rem_list[[1]], but it is not.
## Actual: (bug) it survives in many GenomicRanges versions because of padded vs unpadded f2
out <- GenomicRanges::setdiff(x_list, rem_list, ignore.strand = TRUE)
out[[1]]

## SDeconstruction uses padded f2 based on per-object max(f2) shows the issue
xgr <- GenomicRanges:::deconstructGRLintoGR(x_list)
ygr <- GenomicRanges:::deconstructGRLintoGR(rem_list)

xgr[1]  # "1|04"
ygr[1]  # "1|4"  (note: no zero padding)

The problematic code line I think is in deconstructGRLintoGR:

f12 <- .paste12(f1, f2)