Cannot use findOverlaps with 2 GInteractions objects with different ranges.
1
0
Entering edit mode
Sarah ▴ 10
@c3e29053
Last seen 10 months ago
United States

Hello!

I have a set of genome-wide interactions (gi1). I then subset these interactions to one chromosome and used the anchors to generate new interactions (gi2). My problem now is I want to see which of the new interactions I made were in the original genome-wide set, but since gi1 contains regions from all chromosomes, it has many more ranges than gi2, and R Studio crashes when I try to compare the 2 sets.

Through some testing, I have realized that when 2 GInteractions objects have different ranges in their regions, using either gi1 %in% gi2 or findOverlaps(gi1, gi2) crashes my RStudio session with the "R Session Aborted. R encountered a fatal error." popup.

Here is a simple example that causes the error:

library(InteractionSet)

## Setup modified from the example constructor for GInteractions
# Create first GInteractions object
N <- 30
all.starts <- sample(1000, N)
all.ends <- all.starts + round(runif(N, 5, 20))
all.regions <- GRanges(rep(c("chrA", "chrB"), c(N-10, 10)), IRanges(all.starts, all.ends))

Np <- 20
all.anchor1 <- sample(N, Np)
all.anchor2 <- sample(N, Np)
gi1 <- GInteractions(all.anchor1, all.anchor2, all.regions)

# re-sample to get new ranges
all.starts <- sample(1000, N)
all.ends <- all.starts + round(runif(N, 5, 20))
all.regions <- GRanges(rep(c("chrA", "chrB"), c(N-10, 10)), IRanges(all.starts, all.ends))

all.anchor1 <- sample(N, Np)
all.anchor2 <- sample(N, Np)
gi2 <- GInteractions(all.anchor1, all.anchor2, all.regions)

## Compare interactions
## This line crashes my RStudio session
gi1 %in% gi2

## Alternatively, this line also crashes my RStudio session
findOverlaps(gi1, gi2)

I've run this simple example just using R on the command line and got the following traceback:

 *** caught bus error ***
address 0x100000cfeedfacf, cause 'invalid alignment'

Traceback:
 1: .paired_overlap_finder2(query, subject, select = select, maxgap = maxgap,     minoverlap = minoverlap, type = type, ignore.strand = ignore.strand,     ..., use.region = use.region)
 2: .local(query, subject, maxgap, minoverlap, type, select, ...)
 3: findOverlaps(x, table, type = "equal", select = "first", use.region = "same")
 4: findOverlaps(x, table, type = "equal", select = "first", use.region = "same")
 5: match(x, table, nomatch = 0L)
 6: match(x, table, nomatch = 0L)
 7: gi1 %in% gi2
 8: gi1 %in% gi2

My session info:

> sessionInfo()

R version 4.3.2 (2023-10-31)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Big Sur 11.6

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/New_York
tzcode source: internal

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] InteractionSet_1.30.0       SummarizedExperiment_1.32.0
 [3] Biobase_2.62.0              MatrixGenerics_1.14.0      
 [5] matrixStats_1.2.0           GenomicRanges_1.54.1       
 [7] GenomeInfoDb_1.38.5         IRanges_2.36.0             
 [9] S4Vectors_0.40.2            BiocGenerics_0.48.1        

loaded via a namespace (and not attached):
 [1] SparseArray_1.2.3       zlibbioc_1.48.0         Matrix_1.6-5           
 [4] lattice_0.22-5          abind_1.4-5             GenomeInfoDbData_1.2.11
 [7] S4Arrays_1.2.0          XVector_0.42.0          RCurl_1.98-1.14        
[10] bitops_1.0-7            grid_4.3.2              DelayedArray_0.28.0    
[13] compiler_4.3.2          tools_4.3.2             Rcpp_1.0.12            
[16] crayon_1.5.2

Any help/advice would be appreciated!

InteractionSet GInteractions • 954 views
ADD COMMENT
0
Entering edit mode

More info: This seems to be a problem only on 4.3.2. This code works fine on 4.3.2 with the same package version of InteractionSet. I know there is also a platform difference (mac vs linux), so if anyone else is willing to test the example and let me know if it works for you, that would be helpful too! Thanks!

Session info where the code works:

R version 4.3.1 (2023-06-16)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux 8.8 (Ootpa)

Matrix products: default
BLAS/LAPACK: /nas/longleaf/rhel8/apps/r/4.3.1/lib/libopenblas_zenp-r0.3.23.so;  LAPACK version 3.11.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
 [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: America/New_York
tzcode source: system (glibc)

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] InteractionSet_1.30.0       SummarizedExperiment_1.32.0 Biobase_2.62.0             
 [4] MatrixGenerics_1.14.0       matrixStats_1.1.0           GenomicRanges_1.54.1       
 [7] GenomeInfoDb_1.38.1         IRanges_2.36.0              S4Vectors_0.40.1           
[10] BiocGenerics_0.48.1        

loaded via a namespace (and not attached):
 [1] SparseArray_1.2.2       zlibbioc_1.48.0         Matrix_1.6-3            lattice_0.22-5         
 [5] abind_1.4-5             GenomeInfoDbData_1.2.11 S4Arrays_1.2.0          XVector_0.42.0         
 [9] RCurl_1.98-1.13         bitops_1.0-7            grid_4.3.1              DelayedArray_0.28.0    
[13] compiler_4.3.1          rstudioapi_0.15.0       tools_4.3.1             Rcpp_1.0.11            
[17] crayon_1.5.2
ADD REPLY
0
Entering edit mode
Aaron Lun ★ 28k
@alun
Last seen 10 hours ago
The city by the bay

Hm. Smells like an issue with the C++ code in InteractionSet, though I don't see any issues running it on my own computer (Ubuntu Linux, R-devel). Nothing comes up with valgrind, either.

Perhaps you could debug(InteractionSet:::.paired_overlap_finder2) and confirm exactly where it crashes. The primary suspect would be the .Call to cxx_paired_olaps.

ADD COMMENT
0
Entering edit mode

Thank you so much for the reply! Yes, that's exactly the line that is causing an issue:

*** caught bus error ***
address 0x100000cfeedfacf, cause 'invalid alignment'

Traceback:
 1: .Call(cxx_paired_olaps, left.a1 - 1L, left.a2 - 1L, left.bounds$first -     1L, left.bounds$last, olap$ranges.dex - 1L, right.bounds1$first -     1L, right.bounds1$last, o1 - 1L, right.bounds2$first - 1L,     right.bounds2$last, o2 - 1L, .decode_region_mode(use.region,         c("both", "same", "reverse")), select)
 2: .paired_overlap_finder2(query, subject, select = select, maxgap = maxgap,     minoverlap = minoverlap, type = type, ignore.strand = ignore.strand,     ..., use.region = use.region)
 3: .local(query, subject, maxgap, minoverlap, type, select, ...)
 4: findOverlaps(gi1, gi2)
 5: findOverlaps(gi1, gi2)

Your mention of the C++ code reminded me of a similar issue I was having with the rhdf5 package. I think it's something on my end, because when I installed InteractionSet from source now it's working!

BiocManager::install("InteractionSet", force=TRUE, type = "source")

ADD REPLY

Login before adding your answer.

Traffic: 946 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6