Question: Hello biostar followers,

I use normally "findOverlaps" function in R to select reads of my interest from my BAM files (reading in R with library GenomicAlignments). This times I received a weird results :

> single_read.gal

GAlignments object with 1 alignment and 5 metadata columns:

seqnames strand cigar qwidth start end width njunc | qname seq NH HI nM
|

[1] chrY - 51M 51 42301566 42301616 51 0 | HISEQ:125:CAWV3ACXX:6:1107:18815:33205 TCCAGTCCTG...TGTAAGCCAA 1 1 1
-------
seqinfo: 66 sequences from an unspecified genome

> genomic_region.gr

GRanges object with 1 range and 1 metadata column:

seqnames ranges strand | GeneID
|
[1] chrY 42301297-42301688 + | MTA_Mm
-------
seqinfo: 36 sequences from an unspecified genome; no seqlengths

So, the read "single_read.gal" is clearly within the GRanges interval ( 42301297 < 42301566 & 42301616 < 42301688 ), opposite strand, overlapping width = 51 (width of the read).

But, when I run the "findOverlaps" function ignoring the strand info, there is no overlap fulfilling the input requirements :

> findOverlaps( query = single_read.gal, subject = genomic_region.gr, type = "within", minoverlap = 25, ignore.strand = T )

Hits object with 0 hits and 0 metadata columns:
queryHits subjectHits

-------
queryLength: 1 / subjectLength: 1

It works only with "minoverlap" parameter = 0

> findOverlaps( query = single_read.gal, subject = genomic_region.gr, type = "within", minoverlap = 0, ignore.strand = T )

Hits object with 1 hit and 0 metadata columns:
queryHits subjectHits

[1] 1 1
-------
queryLength: 1 / subjectLength: 1

For any-type overlap works even the former setup :

> findOverlaps( query = single_read.gal, subject = genomic_region.gr, type = "any", minoverlap = 25, ignore.strand = T )

Hits object with 1 hit and 0 metadata columns:
queryHits subjectHits

[1] 1 1
-------
queryLength: 1 / subjectLength: 1

Is that a bug ? Thanks in advance for any respond

> sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 9 (stretch)

Matrix products: default
BLAS: /usr/lib/atlas-base/atlas/libblas.so.3.0
LAPACK: /usr/lib/atlas-base/atlas/liblapack.so.3.0

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

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

other attached packages:
[1] data.table_1.11.4 GenomicAlignments_1.16.0 Rsamtools_1.32.3 Biostrings_2.48.0 XVector_0.20.0 SummarizedExperiment_1.10.1 DelayedArray_0.6.5
[8] BiocParallel_1.14.2 matrixStats_0.54.0 Biobase_2.40.0 GenomicRanges_1.32.6 GenomeInfoDb_1.16.0 IRanges_2.14.11 S4Vectors_0.18.3
[15] BiocGenerics_0.26.0 BiocInstaller_1.30.0

loaded via a namespace (and not attached):
[1] lattice_0.20-35 bitops_1.0-6 grid_3.5.1 zlibbioc_1.26.0 Matrix_1.2-14 tools_3.5.1 RCurl_1.95-4.11 compiler_3.5.1
[9] GenomeInfoDbData_1.1.0

Answer: Seems like a possible bug but also looks to be related to the classes used. If you convert the GAlignments object to GRanges, it works as expected.

## setting up your same objects
genomic_region.gr <- GRanges("chrY:42301297-42301688:+")
single_read.gal <- as(GRanges("chrY:42301566-42301616:-"),'GAlignments')

## test as is
findOverlaps( query = single_read.gal, subject = genomic_region.gr, type = "within", minoverlap = 25, ignore.strand = T )

Hits object with 0 hits and 0 metadata columns:
queryHits subjectHits

-------
queryLength: 1 / subjectLength: 1

## test after changing class
single_read.gal <- as(single_read.gal,'GRanges')

findOverlaps( query = single_read.gal, subject = genomic_region.gr, type = "within", minoverlap = 25, ignore.strand = T )

Hits object with 1 hit and 0 metadata columns:
queryHits subjectHits

[1] 1 1
-------
queryLength: 1 / subjectLength: 1

You might get a better explanation in the Bioconductor forums.