GRanges to subset by overlap within 200kb of start or end of a gene

user_11706046 0
Viewed 100

I have two GRanges datas and I would like to subset them by overlaps such that the overlap could also be present within 200kb of either start or end of the gene.

I was using the following command subsetByOverlaps(gr2, gr1, type = "equal", maxgap = 200000)

using type = "equal" and maxgap= 200000 to get the result I want and I was wondering if it is the correct way to answer my question.

I am not sure if I completely understand the usage of maxgap and hence would like your help or any suggestion in order to get the desired result.

Thanks in advance

Best, S

bioconductor r
1 Answers

I think the proper way to go is subsetByOverlaps(gr2, gr1, type = "any", maxgap = 200000).

From the findOverlaps documentation: "The ‘maxgap’ parameter has special meaning with the special overlap types. For ‘start’, ‘end’, and ‘equal’, it specifies the maximum difference in the starts, ends or both, respectively. For ‘within’, it is the maximum amount by which the subject may be wider than the query. If ‘maxgap’ is set to -1 (the default), it's replaced internally by 0."

In fact, let us consider two IRanges objects as follows:

library(IRanges)

d1 <- IRanges(start = c(1,    2100, 5000, 8000), width = c(100, 300, 400, 600))
d2 <- IRanges(start = c(1000, 2000, 4000, 9000), width = c(300, 200, 400, 100))

d1
##> IRanges object with 4 ranges and 0 metadata columns:
##>           start       end     width
##>       <integer> <integer> <integer>
##>   [1]         1       100       100
##>   [2]      2100      2399       300
##>   [3]      5000      5399       400
##>   [4]      8000      8599       600
d2
##> IRanges object with 4 ranges and 0 metadata columns:
##>          start       end     width
##>      <integer> <integer> <integer>
##>  [1]      1000      1299       300
##>  [2]      2000      2199       200
##>  [3]      4000      4399       400
##>  [4]      9000      9099       100

If we look for intersections between intervals within a gap, with type set either to "equal" or "any":

findOverlapPairs(d1,d2, maxgap = 650, type = "equal")
##>  > Pairs object with 1 pair and 0 metadata columns:
##>           first    second
##>       <IRanges> <IRanges>
##>   [1] 2100-2399 2000-2199


findOverlapPairs(d1,d2, maxgap = 650, type = "any")

##> > Pairs object with 3 pairs and 0 metadata columns:
##>          first    second
##>      <IRanges> <IRanges>
##>  [1] 2100-2399 2000-2199
##>  [2] 5000-5399 4000-4399
##>  [3] 8000-8599 9000-9099

We can see that using "equal" we miss some matches.

user_32075090
Related