Hi,

I know findOverlaps() from GenomicRanges package does pretty good job for finding overlapped regions.I have studied whole vignette of GenomicRanges, BiocParallel pakages. findOverlaps function is very well done, but I need element wise operation across several GRanges objects simultaneously. I come up following reproducible example and put my desire output as well.

Objective: find out overlapped regions across multiple GRanges objects simultanously (a.k.a, element-wise), aim to provide co-localization test and to save discarded enriched regions that was discarded in single trail, but discarded enriched regions will be saved by combining evidence of multiple Chip-seq experiment sample (find out its overlapped regions across across multiple GRanges objects in parallel). 

for example:

Step 1: setup

library(GenomicRanges)
a <- GRanges(
  seqnames=Rle(c("chr1", "chr2", "chr3", "chr4"), c(3, 2, 1, 2)),
  ranges=IRanges(seq(1, by=5, len=8), seq(4, by=5, len=8)),
  rangeName=letters[seq(1:8)],
  score=sample(1:20, 8, replace = FALSE))

b <- GRanges(
  seqnames=Rle(c("chr1", "chr2", "chr3","chr4"), c(4, 1, 2, 1)),
  ranges=IRanges(seq(2, by=6, len=8), seq(5, by=6, len=8)),
  rangeName=letters[seq(1:8)],
  score=sample(1:20, 8, replace = FALSE))

c <- GRanges(
  seqnames=Rle(c("chr1", "chr2", "chr3","chr4"), c(2, 4, 1,1)), 
  ranges=IRanges(seq(4, by=4, len=8), seq(7, by=4, len=8)),
  rangeName=letters[seq(1:8)],
  score=sample(1:15, 8, replace = FALSE))

Step 2: I want to do overlap operation simultaneously such as:

# each iteration, only take one ranges out of multiple ranges in
# querySample (a.k.a, first GRanges objects)

queryRange <- a[1]
ov_b <- b[subjectHits(findOverlaps(a[1], b))]
ov_c <- c[subjectHits(findOverlaps(a[1], c))]

Step 3: then, I want to keep only one overlapped ranges from ov_b, ov_c (if there are multiple intersection were found), such as:

ov_b_keepOne <- ov_b[which.max(ov_b$score)]  # if multiple overlapped
                                             # found, only keep most
                                             # significant one.
ov_c_KeepOne <- ov_c[which.max(ov_c$score)]

Step 4: then, use chisq.test to find out combined pvalue, such as:

comb.pval <- chisq.test(c(queryRange$score, ov_b$score, ov_c$score))$p.value

Step 5: then, given threshold value gamma=1e-8,  (for example, comb.pvalue > gamma), then, execute below:

res <- GRangesList('querySample'=queryRange,
                   'targetSample_1'=ov_b_keepOne,
                   'targetSample_2'=ov_c_KeepOne)

Step 6: finally, this is my expected results :

res.df <- as.data.frame(res)
write.table(res.df, file="Confirmed.bed")

Objective: only take one ranges (a.k.a, each element of GRanges objects) out of multiple ranges, to find its overlapped ranges from multiple GRanges objects simultaneously (seems bplapply function can does this, but I haven't familiar with batch processing in R). I have bottleneck with this problem for a while, still not solved efficiently, any idea or possible approach to save my effort to efficiently develop my packages. Please point me out what should I do? Thanks a lot

Best regards:

Jurat

I updated your question to make the formatting clearer, and to include steps so that we can refer to the individual steps involved.

From the discussion below, an efficient starting point for this work flow seems to be

## helper function: identify the overlapping range satisfying FUN
keepone <- function(gr, hitlist, FUN=which.max) {
    idx0 <- as(FUN(extractList(gr$score, hitlist)), "List")
    idx1 <- unlist(extractList(seq_along(gr), hitlist)[idx0])
    ## FIXME: what about NA's when there are no matching ranges?
    gr[idx1]
}

ab_hitlist <- as(findOverlaps(a, b), "List")
ac_hitlist <- as(findOverlaps(a, c), "List")

test <- (lengths(ab_hitlist) > 0) + (lengths(ac_hitlist) > 0)
keep <- test >= minNumOlp_Required
a = a[keep]
ab_hitlist = ab_hitlist[keep]
ac_hitlist = ac_hitlist[keep]

b_one <- keepone(b, ab_hitlist)
c_one <- keepone(c, ac_hitlist)

comb.pval <-
    mapply(function(a, b, c) chisq.test(c(a, b, c))$p.value,
           a$score, b_one$score, c_one$score)

Updating each step

The place I want to start is with step 2, which I will replace with

ab_hitlist <- as(findOverlaps(a, b), "List")
ac_hitlist <- as(findOverlaps(a, c), "List")

for the reasons below.

I will replace step 4 with

comb.pval  <- mapply(function(a, b, c) {
    chisq.test(c(a, b, c))$p.value
}, a$score, max(extractList(b$score, ab_hitlist)), max(extractList(c$score, ac_hitlist)))

Step 3 is to identify the overlapping range with maximum score. I did this as follows:

keepone <- function(gr, hitlist) {
    idx0 <- as(which.max(extractList(gr$score, hitlist)), "List")
    idx1 <- unlist(extractList(seq_along(gr), hitlist)[idx0])
    ## FIXME: what about NA's when there are no matching ranges?
    gr[idx1]
}
b_one <- keepone(b, ab_hitlist)
c_one <- keepone(c, ac_hitlist)

But now that the specification for step 4 has been clarified as the maximum score of the overlapping ranges, and we have the overlapping range, step 4 is

comb.pval  <- mapply(function(a, b, c) chisq.test(c(a, b, c))$p.value,
                     a$score, b_one$score, c_one$score) 

Justification for updated step 2:

It is very inefficient to calculate overlaps between 1 range and many ranges. So 'invert' the problem and find all overlaps

ab_hit <- findOverlaps(a, b)
ac_hit <- findOverlaps(a, c)

For convenience, split the hits into a list, where each element of the list corresponds to the overlaps for a different query

ab_hitlist <- as(ab_hit, "List")

For your example, this is

> ab_hitlist
IntegerList of length 8
[[1]] 1 2 3 4
[[2]] 1 2 3 4
[[3]] 1 2 3 4
[[4]] 5
[[5]] 5
[[6]] 6 7
[[7]] 8
[[8]] 8

Convince yourself that ab_hitlist[[1]]

> ab_hitlist[[1]]
[1] 1 2 3 4

is the same as subjectHits(findOverlaps(a[1], b)), and so on for a[2], a[3], a[4]. The reason for this change is that taking ranges one-at-a-time scales linearly

> library(microbenchmark)
> f <- function(a, b)
+     lapply(seq_along(a), function(i) subjectHits(findOverlaps(a[i], b)))
> g <- function(a, b)
+     as(findOverlaps(a, b), "List")
> a0 <- a[1:4]; a2 <- c(a, a); a4 <- c(a2, a2); a8 <- c(a4, a4); a16 <- c(a8, a8)
> microbenchmark(f(a0, b), f(a, b), f(a2, b), times=10)
Unit: milliseconds
     expr       min        lq      mean    median        uq       max neval cld
 f(a0, b)  75.21582  76.07495  77.63697  77.76686  78.44264  80.56188    10 a  
  f(a, b) 152.56619 154.25827 157.34376 156.14110 158.94034 169.31483    10  b 
 f(a2, b) 309.08562 309.42790 313.71373 312.01247 316.88869 323.32852    10   c
> microbenchmark(g(a0, b), g(a, b), g(a2, b), g(a4, b), g(a8, b), g(a16, b),
+                times=10)
Unit: milliseconds
      expr      min       lq     mean   median       uq       max neval cld
  g(a0, b) 17.24473 17.32625 17.93150 17.67085 18.01294  20.22646    10   a
   g(a, b) 17.33038 17.50243 18.09587 17.73340 18.19032  20.01569    10   a
  g(a2, b) 20.46892 20.52388 21.30909 21.11872 21.66259  23.64361    10   a
  g(a4, b) 20.51881 20.68018 21.32741 21.17813 21.65589  23.31553    10   a
  g(a8, b) 20.55628 20.64502 21.14069 20.86133 21.33573  22.41963    10   a
 g(a16, b) 20.50725 20.82241 34.97887 21.24104 21.45764 158.24840    10   a

whereas the same is not true for creating a hits object (the scaling is not evident from the small example data).

Justification for updated step 4.

The extractList() function will take a vector-like object as a first argument, and reshape it according to the indexes of the second argument. It is relatively amazingly useful to see that

> extractList(b$score, ab_hitlist)
IntegerList of length 8
[[1]] 12 19 11 6
[[2]] 12 19 11 6
[[3]] 12 19 11 6
[[4]] 7
[[5]] 7
[[6]] 10 20
[[7]] 14
[[8]] 14

extracts the scores and reshapes them into the correct geometry of the hits.

The calculation of p-values applies chisq.test() to the element-wise 'parallel' vectors a$score, and from the discussion in the comments below to max(extractList(b$score, ab_hitlist)), and max(extractList(c$score, ac_hitlist)). So

comb.pval <- mapply(function(a, b, c) {
    chisq.test(c(a, b, c))$p.value
}, a$score, max(extractList(b$score, ab_hitlist)), max(extractList(c$score, ac_hitlist)))

The result is

> comb.pval
[1] 1.380288e-04 1.290620e-04 3.281973e-05 5.494160e-01 4.821936e-01
[6] 1.813352e-01 3.805041e-03 8.590223e-02

All other operations in can be 'vectorized' to be efficient, but chisq.test() cannot be (I bet there's a function that does do this in a vectorized way?). It is therefore a candidate for parallel evaluation, but I would only do this at the end, when it becomes apparent that this is a rate-limiting step.

library(BiocParallel)
comb.pval <- bpmapply(function(a, b, c) {
    chisq.test(c(a, b, c))$p.value
}, a$score, max(extractList(b$score, ab_hitlist)), max(extractList(c$score, ac_hitlist)))

Justification for updated step 3.

which.max(extractList(b$score, ab_hitlist))

returns a plain vector of indexes in to each of the list elements of the hitlist. Coercing this to a List

idx0 <- as(which.max(extractList(gr$score, hitlist)), "List")

restores the 'geometry' of the vector, so that we have a list of length(ab_hitline). We then create an index into the grang that we're interested in, cast it in the appropriate geometry, and subset with the index

idx1 <- unlist(extractList(seq_along(b), ab_hitlist)[idx0])

These are the indexes into the original ranges, so we can retrieve them with

b[idx1]

A little complicated, so we implement this as a small helper function

keepone <- function(gr, hitlist) {
    idx0 <- as(which.max(extractList(gr$score, hitlist)), "List")
    idx1 <- unlist(extractList(seq_along(gr), hitlist)[idx0])
    gr[idx1]
}

I think you want to look into the BiocParallel package. There is also GenomicFiles, which has specific functionality for processing files in parallel.

It would look something like this:

hitsList <- bplapply(grs[-1L], function(gr) findOverlaps(grs[[1L]], gr))

If you need to do things pairwise, look into expand.grid() and bpmapply(). Once you compute the hits in parallel, and you have a list of Hits, do this to intersect them:

hits <- Reduce(intersect, hitsList)

Did you tried ChIPpeakAnno::findOverlapsOfPeaks like this?

library(ChIPpeakAnno)
ol <- findOverlapsOfPeaks(a, b, c, maxgap=0)
olpeaks <- ol$peaklist[["a///b///c"]]
all <- GRangesList(a, b, c)
all <- mapply(function(.ele, .name){names(.ele)<- paste(.name, 1:length(.ele), sep="__"); .ele}, all, c("a", "b", "c"))
all <- unlist(GRangesList(all))
pval <- sapply(olpeaks$peakNames, function(.ele) chisq.test(all[.ele]$score)$p.value) ##here should be changed to keep the best one for each peak list
...
​

You could also have a look and try regioneR and LOLA bioconductor packages.