Hi,

I have a GRanges object und would like to retrieve sets of reciprocal overlaps:

library(GenomicRanges)

gr = GRanges(seqnames = rep("chr1",4),
             ranges   = IRanges(start = c(1,2,3,8),
                                 end  = c(12,6,6,10)),
             strand = rep("+", 4))
> gr
GRanges object with 4 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr1   [1, 12]      +
  [2]     chr1   [2,  6]      +
  [3]     chr1   [3,  6]      +
  [4]     chr1   [8, 10]      +
  -------

Range #1 overlaps with all other ranges, however I would like to get only those overlaps which hold in a transitive manner. such as: Range #1 overlaps with #2, #2 with #3 and #3 with #1.
At the moment I am filtering the output of findOverlaps (unsatisfactorily with the help of a loop..). Is there a more direct way to do this?

Many thanks for any comments,

Stefanie

The problem could be thought of in graph terms.  The nodes are the regions.  The edges consist of connecting all regions that overlap.  These edges are exactly what is returned by findOverlaps().  If you take these concepts and actually construct a graph, there are a number of packages that will construct the transitive closure.  This eliminates the need to do findOverlaps() in the loop.  

Edit: Vince points out that the RBGL package has a transitive closure functionality; this should be pretty performant, I suspect.

Edit: Here is some code with a largish set of regions (10,000 exons) showing performance (about 5 seconds for the entire script, including transitive closure calculation).

# Form a large set of regions
library("TxDb.Hsapiens.UCSC.hg19.knownGene")
exons = exons(TxDb.Hsapiens.UCSC.hg19.knownGene)[1:10000]
library(graph)
# Make a graph with edges representing overlapping regions
g = graphNEL(nodes=as.character(seq_along(exons)))
x = as.data.frame(findOverlaps(exons))
g = addEdge(from=as.character(x[,1]),to=as.character(x[,2]),g)
# Compute the transitive closure
library(RGBL)
system.time((g1 = transitive.closure(g)))
  user  system elapsed 
  2.666   0.426   3.110 
g1
# And use edges(g1) to capture the transitive relationships.
head(edges(g1))
$`1`
[1] "1"

$`2`
[1] "2" "3" "4"

$`3`
[1] "2" "3" "4"

$`4`
[1] "2" "3" "4"

$`5`
[1] "5" "6"

$`6`
[1] "5" "6"

Hi Herve,

I think one way to achieve what I would like to do is as follow:

ov = as.matrix(findOverlaps(gr, gr))
ranges_by_overlap = split(gr[ov[, 2]], ov[, 1])

lapply(ranges_by_overlap, function(x) Reduce(intersect,x))

By doing so, all 'correct' transitive overlaps will have a valid intersection as output from Reduce(intersect,x). All non-transitive overlaps will return an empty range.

However, this is unfortunately pretty slow...
Best,

Stefanie

Hi Stephanie,

Not sure what you're trying to achieve exactly but note that the easiest way to "cluster ranges by overlap" is by calling reduce() with with.revmap=TRUE:

library(GenomicRanges)
gr <- GRanges(seqnames="chr1",
              ranges=IRanges(c( 1, 2, 3,  8, 15, 20,  9),
                             c(12, 6, 6, 10, 20, 25, 13)),
              strand="+")
rgr <- reduce(gr, min.gapwidth=0, with.revmap=TRUE)

Note that min.gapwidth is 1 by default and that means that adjacent ranges that do not overlap would be considered to be in the same cluster. So we use min.gapwidth=0 to avoid that.

Then:

> rgr
GRanges object with 2 ranges and 1 metadata column:
      seqnames    ranges strand |        revmap
         <Rle> <IRanges>  <Rle> | <IntegerList>
  [1]     chr1  [ 1, 13]      + |     1,2,3,...
  [2]     chr1  [15, 25]      + |           5,6
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

rgr contains 1 range per cluster and its revmap metadata column contains the indices of the original ranges that belong to each cluster:

> mcols(rgr)$revmap
IntegerList of length 2
[[1]] 1 2 3 4 7
[[2]] 5 6

Hope this helps,

H.