Answers.Body: <p>If preferred, this can be achieved with base R functions:</p>
  <pre><code># First store the chromosome and start position of your gene of interest in variables
  &gt; chrPos &lt;- as.character(seqnames(gr))[gr<span class="math-container">$GeneSymbol=="gene3"]
  &gt; startPos &lt;- start(ranges(gr))[gr$</span>GeneSymbol==&quot;gene3&quot;]
  
  # 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
  &gt; gr$distance &lt;- 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)
  &gt; grOut &lt;- gr[seqnames(gr) %in% chrPos &amp; (abs(gr$distance) &lt;= 1000000)]
  &gt; grOut
  GRanges object with 3 ranges and 4 metadata columns:
      seqnames         ranges strand |     score        GC  GeneSymbol  distance
         &lt;Rle&gt;      &lt;IRanges&gt;  &lt;Rle&gt; | &lt;integer&gt; &lt;numeric&gt; &lt;character&gt; &lt;integer&gt;
    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
  </code></pre>
  <p>This strategy provides the added benefit of adding a 'distance-from-gene' (<code>gr$distance</code>) column to the dataset which can be useful for visualizations.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-03-04T16:44:37.770Z
Answers.LastActivityDate: 2022-03-04T18:25:12.067Z
Answers.LastEditDate: 2022-03-04T18:25:12.067Z
Answers.UserDisplayName: Gawain
Answers.Body: <p>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:</p>
  <pre><code>library(&quot;GenomicRanges&quot;)
  
  &gt; chrPos &lt;- as.character(seqnames(gr))[gr<span class="math-container">$GeneSymbol=="gene3"]
  &gt; startPos &lt;- start(ranges(gr))[gr$</span>GeneSymbol==&quot;gene3&quot;]
  &gt; range &lt;- 1000000
  
  &gt; q &lt;- GRanges(seqnames=chrPos,
                 ranges=IRanges(start=startPos-range, end=startPos+range))
  &gt; grOut &lt;- subsetByOverlaps(gr, q)
  &gt; grOut
  GRanges object with 3 ranges and 3 metadata columns:
      seqnames         ranges strand |     score        GC  GeneSymbol
         &lt;Rle&gt;      &lt;IRanges&gt;  &lt;Rle&gt; | &lt;integer&gt; &lt;numeric&gt; &lt;character&gt;
    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
  </code></pre>
  <p>Credit to <a href="https://www.biostars.org/p/351377/#351383" rel="nofollow noreferrer">@zx8754</a> for inspiring this solution.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-02-24T01:30:12.287Z
Answers.LastActivityDate: 2022-03-04T18:31:54.143Z
Answers.LastEditDate: 2022-03-04T18:31:54.143Z
Answers.UserDisplayName: Gawain
Body: <p><strong>Goal:</strong> 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.</p>
  <p>I am aware of some functions that come close:</p>
  <ul>
  <li><code>GenomicRanges::nearestKNeighbors()</code> (and related <code>distanceToNearest()</code>, <code>nearest()</code>,
  <code>precede()</code>, <code>follow()</code>)</li>
  <li><code>GenomicRanges::subsetByOverlaps()</code></li>
  <li><code>GenomicRanges::window()</code></li>
  </ul>
  <p>And I am aware of the UNIX command-line 'Bedtools' tools: <code>window</code>, <code>closest</code>, and <code>intersect</code>, however I wish to accomplish this task using R only and the above functions have not been doing exactly as I would like.</p>
  <p>Here is an example GRanges dataset:</p>
  <pre><code>&gt; gr &lt;- GRanges(
          seqnames = Rle(c(&quot;chr1&quot;, &quot;chr2&quot;, &quot;chr1&quot;, &quot;chr3&quot;), 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(&quot;-&quot;, &quot;+&quot;, &quot;*&quot;, &quot;+&quot;, &quot;-&quot;)), c(1, 2, 2, 3, 2)),
          score = 1:10, GC = seq(1, 0, length=10), GeneSymbol = c(paste(&quot;gene&quot;, seq(1:10), sep=&quot;&quot;)))
  
  &gt; gr
  GRanges object with 3 ranges and 3 metadata columns:
      seqnames          ranges strand |     score        GC  GeneSymbol
         &lt;Rle&gt;       &lt;IRanges&gt;  &lt;Rle&gt; | &lt;integer&gt; &lt;numeric&gt; &lt;character&gt;
    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
  </code></pre>
  <p>Ideally, by specifying the name (e.g., <code>&quot;gene3&quot;</code>) or start coordinate (e.g., <code>chr2</code> <em>AND</em> <code>677777</code>) of the gene of interest AND a window size (e.g., <code>1000000</code> or <code>1 Mb</code> upstream <em>AND</em> downstream) to a function, I hope to get an output like below:</p>
  <pre><code>GRanges object with 3 ranges and 3 metadata columns:
      seqnames          ranges strand |     score        GC  GeneSymbol
         &lt;Rle&gt;       &lt;IRanges&gt;  &lt;Rle&gt; | &lt;integer&gt; &lt;numeric&gt; &lt;character&gt;
    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
  </code></pre>
  <p>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 <em>more than 1 Mb away</em> from the start coordinate of gene3.</p>
  <hr />
  <p>I have tried various methods to isolate the rows within 1 Mb of my gene of interest with no luck:</p>
  <pre><code>&gt; gr1 &lt;- gr[grep(&quot;gene3&quot;, gr$GeneSymbol),]
  &gt; grOut &lt;- subsetByOverlaps(gr, gr1, maxgap = 1000000, invert=T)
  
  &gt; grOut &lt;- gr[mcols(gr)$GeneSymbol == &quot;gene3&quot; &amp; start(gr) &lt; 677777-1000000 &amp; end(gr) &gt; 1877777+1000000]
  
  &gt; grOut &lt;- window(gr, GRanges(&quot;chr2&quot;, IRanges(start=677777-1000000, end=1877777+1000000)))
  Error: 'start' must be a single number or NA
  
  &gt; grOut &lt;- gr[IRanges(start=677777, end=1877777)]
  Error: subscript contains out-of-bounds ranges
  
  &gt; grOut &lt;- nearest(gr, GRanges(&quot;chr2&quot;, IRanges(start=677777-1000000, end=1877777+1000000)), select=&quot;all&quot;, ignore.strand=T)
  </code></pre>
  <p>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?</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-03-04T18:52:05.980Z
Comments.Text: For those wishing to find/verify the coordinates for a gene of interest (GOI) in their GRanges object **when only the gene symbol/ID/name is known**: use `gr[grep("<GOI>", names(gr)),]` when *column names* have been set as gene names ... use `gr[grep("<GOI>", rownames(gr)),]` when *row names* have been set as gene names ... use `gr[grep("<GOI>", mcols(gr)$<geneNameColumn>),]` when a *metadata column* contains gene names.
Comments.UserDisplayName: Gawain
Tags: <r><filtering><genomicranges>
Title: How to subset a GRanges object based on a specific genomic window of interest?
ThreadText: # How to subset a GRanges object based on a specific genomic window of interest?
  
  (Asked by: username_0 on 2022-02-23)
  
  **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?
  
  **Tags:** r, filtering, genomicranges
  
  > **Comments:**
  > **username_0**: For those wishing to find/verify the coordinates for a gene of interest (GOI) in their GRanges object **when only the gene symbol/ID/name is known**: use `gr[grep("<GOI>", names(gr)),]` when *column names* have been set as gene names ... use `gr[grep("<GOI>", rownames(gr)),]` when *row names* have been set as gene names ... use `gr[grep("<GOI>", mcols(gr)$<geneNameColumn>),]` when a *metadata column* contains gene names.  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_0 on 2022-02-24. Score: 4 | Accepted answer
  
  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](https://www.biostars.org/p/351377/#351383) for inspiring this solution.
  
  
  ---
  
  ### Answer by username_0 on 2022-03-04. Score: 0
  
  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.
  
  
  ---
  
  