How to subset a GRanges object based on a specific genomic window of interest?

Question

Goal: To keep only rows of a GRanges object if their start coordinate is within 1 Mb of the start coordinate of a gene of interest.

I am aware of some functions that come close:

• GenomicRanges::nearestKNeighbors() (and related distanceToNearest(), nearest(), precede(), follow())
• GenomicRanges::subsetByOverlaps()
• GenomicRanges::window()

And I am aware of the UNIX command-line 'Bedtools' tools: window, closest, and intersect, however I wish to accomplish this task using R only and the above functions have not been doing exactly as I would like.

Here is an example GRanges dataset:

> gr <- GRanges(
        seqnames = Rle(c("chr1", "chr2", "chr1", "chr3"), c(1, 3, 2, 4)),
        ranges = IRanges(start = seq(300000, 2000000, length=10), end = seq(700000, 6000000, length=10), names = head(letters, 10)),
        strand = Rle(strand(c("-", "+", "*", "+", "-")), c(1, 2, 2, 3, 2)),
        score = 1:10, GC = seq(1, 0, length=10), GeneSymbol = c(paste("gene", seq(1:10), sep="")))

> gr
GRanges object with 3 ranges and 3 metadata columns:
    seqnames          ranges strand |     score        GC  GeneSymbol
       <Rle>       <IRanges>  <Rle> | <integer> <numeric> <character>
  a     chr1   300000-700000      - |         1  1.000000       gene1
  b     chr2  488888-1288888      + |         2  0.888889       gene2
  c     chr2  677777-1877777      + |         3  0.777778       gene3
  d     chr2  866666-2466666      * |         4  0.666667       gene4
  e     chr1 1055555-3055555      * |         5  0.555556       gene5
  f     chr1 1244444-3644444      + |         6  0.444444       gene6
  g     chr3 1433333-4233333      + |         7  0.333333       gene7
  h     chr3 1622222-4822222      + |         8  0.222222       gene8
  i     chr3 1811111-5411111      - |         9  0.111111       gene9
  j     chr3 2000000-6000000      - |        10  0.000000      gene10
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

Ideally, by specifying the name (e.g., "gene3") or start coordinate (e.g., chr2 AND 677777) of the gene of interest AND a window size (e.g., 1000000 or 1 Mb upstream AND downstream) to a function, I hope to get an output like below:

GRanges object with 3 ranges and 3 metadata columns:
    seqnames          ranges strand |     score        GC  GeneSymbol
       <Rle>       <IRanges>  <Rle> | <integer> <numeric> <character>
  b     chr2  488888-1288888      + |         2  0.888889       gene2
  c     chr2  677777-1877777      + |         3  0.777778       gene3
  d     chr2  866666-2466666      * |         4  0.666667       gene4
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

Notice that only three rows from the original GRanges object remain in this example output because the start coordinates for all other genes were either (i) on different chromosomes or (ii) were more than 1 Mb away from the start coordinate of gene3.

---

I have tried various methods to isolate the rows within 1 Mb of my gene of interest with no luck:

> gr1 <- gr[grep("gene3", gr$GeneSymbol),]
> grOut <- subsetByOverlaps(gr, gr1, maxgap = 1000000, invert=T)

> grOut <- gr[mcols(gr)$GeneSymbol == "gene3" & start(gr) < 677777-1000000 & end(gr) > 1877777+1000000]

> grOut <- window(gr, GRanges("chr2", IRanges(start=677777-1000000, end=1877777+1000000)))
Error: 'start' must be a single number or NA

> grOut <- gr[IRanges(start=677777, end=1877777)]
Error: subscript contains out-of-bounds ranges

> grOut <- nearest(gr, GRanges("chr2", IRanges(start=677777-1000000, end=1877777+1000000)), select="all", ignore.strand=T)

The last attempt recorded above was likely the closest I've gotten to my desired output by storing the index of rows in the GRanges object falling within the range specified, but it presents this index in the form of a SortedByQueryHits object, not a GRanges object as I require. It seems simple enough a task, but I haven't solved it yet; any thoughts or solutions?

Answer 1

Specifying the range of interest (+/- 1 Mb of start coordinate) as a new GRanges object and then subsetting the original GRanges object based on the overlap of genes with the query works:

library("GenomicRanges")

> chrPos <- as.character(seqnames(gr))[gr$GeneSymbol=="gene3"]
> startPos <- start(ranges(gr))[gr$GeneSymbol=="gene3"]
> range <- 1000000

> q <- GRanges(seqnames=chrPos,
               ranges=IRanges(start=startPos-range, end=startPos+range))
> grOut <- subsetByOverlaps(gr, q)
> grOut
GRanges object with 3 ranges and 3 metadata columns:
    seqnames         ranges strand |     score        GC  GeneSymbol
       <Rle>      <IRanges>  <Rle> | <integer> <numeric> <character>
  b     chr2 488888-1288888      + |         2  0.888889       gene2
  c     chr2 677777-1877777      + |         3  0.777778       gene3
  d     chr2 866666-2466666      * |         4  0.666667       gene4
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

Credit to @zx8754 for inspiring this solution.

Answer 2

If preferred, this can be achieved with base R functions:

# First store the chromosome and start position of your gene of interest in variables
> chrPos <- as.character(seqnames(gr))[gr$GeneSymbol=="gene3"]
> startPos <- start(ranges(gr))[gr$GeneSymbol=="gene3"]

# Then add a column to the GRanges object with distances from the start position of each
# gene in the GRanges object to the start position of your gene of interest
> gr$distance <- start(ranges(gr))-startPos

# Finally subset the GRanges object based on your desired window of interest (must be
# symmetrical —— e.g., +/- 1 Mb —— for the following to work)
> grOut <- gr[seqnames(gr) %in% chrPos & (abs(gr$distance) <= 1000000)]
> grOut
GRanges object with 3 ranges and 4 metadata columns:
    seqnames         ranges strand |     score        GC  GeneSymbol  distance
       <Rle>      <IRanges>  <Rle> | <integer> <numeric> <character> <integer>
  b     chr2 488888-1288888      + |         2  0.888889       gene2   -188889
  c     chr2 677777-1877777      + |         3  0.777778       gene3         0
  d     chr2 866666-2466666      * |         4  0.666667       gene4    188889
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

This strategy provides the added benefit of adding a 'distance-from-gene' (gr$distance) column to the dataset which can be useful for visualizations.