I have 2 GRanges objects of equal length.  For each individual pair of rows from the 2 objects, I want to get:

    a) the intersection of the ranges

    b) the ranges unique to the first range

    c) the ranges unique the second range

See my code with comments for specifics, here I will present an overview of my problem.  I have no problem with intersect() - the seems to work.  My trouble is getting the asymmetric differences between paired ranges.  I can do this one row at a time using setdiff() on same row from both elements, but the runtime explodes.  If I pass both GRanges objects, ranges that happen to be overlapping are collapsed into 1 range, which is not what I want.  Because some ranges contain their counterpart from the other GRanges object, psetdiff() crashes.  What is the fastest way to get what I want?

# test data

 > test.a <- GRanges(seqnames=rep("chr1", times=5), ranges=IRanges(start=c(1000, 2000, 4000, 6000, 7000), end=c(2000, 3000, 5000, 7000, 8000)))
 > test.b <- GRanges(seqnames=rep("chr1", times=5), ranges=IRanges(start=c(1000, 2100, 4500, 5900, 7700), end=c(1500, 2900, 5000, 6300, 8200)))

#  the output I am looking for

# note that the results for 1 and 2 share a base in common

# I still count these as 2 separate results

 setdiff(test.a[1], test.b[1])
 GRanges object with 1 range and 0 metadata columns:
       seqnames       ranges strand
          <Rle>    <IRanges>  <Rle>
   [1]     chr1 [1501, 2000]      *
   -------
   seqinfo: 1 sequence from an unspecified genome; no seqlengths

setdiff(test.a[2], test.b[2])
 GRanges object with 2 ranges and 0 metadata columns:
       seqnames       ranges strand
          <Rle>    <IRanges>  <Rle>
   [1]     chr1 [2000, 2099]      *
   [2]     chr1 [2901, 3000]      *
   -------
   seqinfo: 1 sequence from an unspecified genome; no seqlengths



# when i pass the full GRanges objects

# note that 1 and 2 are now combined

# this is not what i want

setdiff(test.a, test.b, resolve.empty="none")
 GRanges object with 4 ranges and 0 metadata columns:
       seqnames       ranges strand
          <Rle>    <IRanges>  <Rle>
   [1]     chr1 [1501, 2099]      *
   [2]     chr1 [2901, 3000]      *
   [3]     chr1 [4000, 4499]      *
   [4]     chr1 [6301, 7699]      *
   -------
   seqinfo: 1 sequence from an unspecified genome; no seqlengths



# calling psetdiff(), which i understand is supposed to be the element-wise setdiff()

# it crashes

# this is not what i want

psetdiff(test.a, test.b, resolve.empty="none")
 Error in start(value) :
   error in evaluating the argument 'x' in selecting a method for function 'start': Error in psetdiff(extractROWS(ranges(x), ok), extractROWS(ranges(y), ok)) :
   some ranges in 'y' have their end points strictly inside
   the range in 'x' that they need to be subtracted from.
   Cannot subtract them.

The psetdif() operation is indeed what you want; however, it is what we call an "endomorphism", i.e., it returns an object of the same type as its inputs. A GRanges is obviously incapable of representing the result of a setdiff when "x" contains "y".  Thus, we need a data structure that is capable of representing that complexity. That data structure is GRangesList.

The solution is to convert your GRanges objects to GRangesList objects and then call psetdiff().

grl.a <- as(test.a, "List")
grl.b <- as(test.b, "List")
psetdiff(grl.a, grl.b)

Hope that helps.