Answers.Body: <p>Here's a way to use <a href="https://github.com/bedops/bedops" rel="nofollow noreferrer">BEDOPS</a>, which was designed to work fast by using sorted input. Other tools now use sorting to accomplish similar performance benefits.</p>
  
  <ol>
  <li><p>Convert GTF annotations to a sorted BED file of genes:</p>
  
  <pre><code><span class="math-container">$ awk '($</span>3=="gene")' annotations.gtf | gtf2bed - &gt; genes.bed
  </code></pre></li>
  <li><p>Sort your intervals, if unsorted:</p>
  
  <pre><code>$ sort-bed intervals.unsorted.bed &gt; intervals.bed
  </code></pre></li>
  <li><p>Filter intervals against 1kb-padded genes:</p>
  
  <pre><code>$ bedops --merge genes.bed | bedops --range 1000 --everything - | bedops --not-element-of 1 intervals.bed - &gt; intervals.filtered.bed
  </code></pre></li>
  </ol>
  
  <p>We first <code>--merge</code> genes to simplify dealing with overlapping genes and further speed downstream filtering.</p>
  
  <p>The file <code>intervals.filtered.bed</code> will contain intervals from <code>intervals.bed</code> which lie outside of a 1kb-window around merged gene space. </p>
  
  <hr>
  
  <p>To run all of this as a one-liner (assuming use of the <code>bash</code> shell):</p>
  
  <pre><code><span class="math-container">$ bedops --merge &lt;(awk '($</span>3=="gene")' annotations.gtf | gtf2bed -) | bedops --range 1000 --everything - | bedops --not-element-of 1 &lt;(sort-bed intervals.unsorted.bed) - &gt; intervals.filtered.bed
  </code></pre>
  
  <p>This creates the same output result, but avoids creation of intermediate files and so should run even faster.</p>
  
  <hr>
  
  <p><strong>Note</strong>: BEDOPS is agnostic about chromosome sizes — this toolkit works with any genome assembly. However, using <code>--range</code> takes care of ranges that would have a negative start position, setting their start coordinates to <code>0</code>. Bounds aren't really relevant here, in any case, because the <code>bedops --not-element-of</code> step queries input intervals, which should already be within chromosome bounds, unless they were generated incorrectly.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-11-01T19:16:54.723Z
Answers.Text: I think there's a missing quote on the first `awk` call
Answers.UserDisplayName: gc5
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-11-02T21:40:23.507Z
Answers.Text: When executing 1 I get `Error: Potentially missing gene or transcript ID from GTF attributes (malformed GTF at line [1]?)`
Answers.UserDisplayName: gc5
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-11-02T23:01:28.997Z
Answers.Text: https://www.biostars.org/p/206342/#206350
Answers.UserDisplayName: Alex Reynolds
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-10-29T17:08:16.627Z
Answers.LastActivityDate: 2018-11-02T14:42:35.780Z
Answers.LastEditDate: 2018-11-02T14:42:35.780Z
Answers.UserDisplayName: Alex Reynolds
Answers.Body: <p>In one line, using <code>bedtools</code></p>
  
  <pre><code>zcat  Homo_sapiens.GRCh38.93.gtf.gz \
   | awk '$3=="gene"' \
   | bedtools slop -b 10000 -g contigs.tsv -i - \
   | bedtools intersect -u -a intervals.bed -b - 
  </code></pre>
  
  <p>This first takes the genes and filters them for the third column being gene. This is important because all <code>GTF</code> lines contain the word <code>gene</code> as <code>gene_id</code> is a manditory attribute for <code>GTF</code>.</p>
  
  <p>Next it passes it to slop, which adds 1 kb either side (ensuring not to go below 0 or off the end of the chr: this isn't that important because we are never outputting these regions anyway, but <code>bedtools</code> requires it). We tell <code>bedtools</code> to work on the standard in by using <code>-i -</code>.</p>
  
  <p>Finally we intersect these gene regions with the <code>intervals.bed</code>. <code>-u</code> says output the line from <code>-a</code> if it overlaps with at least one entry in <code>-b</code> and passing <code>-b -</code> means read file b from the standard in.</p>
  
  <p>This will only work if your <code>GTF</code> contains <code>gene</code> lines, which the GENCODE/Ensembl ones do. However, it is interesting to note that these GTFs are not actaully valid GTF files, as to be valid a GTF file must have both a <code>gene_id</code> and <code>transcript_id</code> attribute, and the gene lines in the GENCODE GTFs to not have this. </p>
  
  <p>If your GTF does not contain gene lines, but can create them by merging all the exons in all the transcripts and filling in the gaps using <code>gtf2gtf</code> from <code>cgat</code>.</p>
  
  <pre><code> cgat gtf2gtf --method=merge-transcripts --stdin Homo_sapiens.GRCh38.93.gtf.gz \
   | bedtools slop -b 10000 -g contigs.tsv -i - \
   | bedtools intersect -u -a intervals.bed -b -
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-11-01T19:58:53.137Z
Answers.LastActivityDate: 2018-11-04T12:51:52.710Z
Answers.LastEditDate: 2018-11-04T12:51:52.710Z
Answers.UserDisplayName: Ian Sudbery
Answers.Body: <p>First I prepared a bed file in which the gene intervals are augmented by 1KB before and after the gene start and end coordinates. Then I intersected this bed file with my original one with the option <code>-wa</code>, therefore retaining only the intervals in my original bed file that intersect with the bedfile produced from the original gtf, filtered by gene regions. The <code>-u</code> flag is also used to print only one line per overlap between the two files, in case an interval in the first overlaps with multiple intervals in the second.</p>
  
  <pre><code># download information about chromosomes
  wget http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.chrom.sizes ./
  
  # Original GTF has a different notation than UCSC
  cat Homo_sapiens.GRCh38.93.gtf | \
    sed -e <span class="math-container">$'s/^./chr&amp;/' -e $</span>'s/^chrMT/chrM/' | \
    awk '<span class="math-container">$3=="gene"' | \
    bedtools slop -b 1000 -g hg38.chrom.sizes -i - | \
    sed -e $</span>'s/\tgene_id.*gene_name/\t/' -e <span class="math-container">$'s/gene_source.*//' | \
    tr -d '"; ' | \
    awk -v OFS=$</span>'\t' '{print <span class="math-container">$1, $</span>4, <span class="math-container">$5, $</span>9}' &gt; genes_1kb.bed
  
  bedtools intersect -wa -u -a intervals.bed -b genes_1kb.bed &gt; intervals.1kb.bed
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-10-29T12:26:00.107Z
Answers.Text: note that this will filter out intervals that overlap a gene, or 1kp upstream, but will leave features that are within 1kbp downstream. And it will only work if the GTF has gene lines, which many GTF files don't have.
Answers.UserDisplayName: Ian Sudbery
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-10-30T21:18:34.110Z
Answers.Text: @IanSudbery, why this will filter out intervals that overlap a gne or 1kb upstream?
Answers.UserDisplayName: gc5
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-10-30T22:38:38.127Z
Answers.Text: Oops I noticed now I forgot to add 1kb to the end. Thanks! Does it solve the problem?
Answers.UserDisplayName: gc5
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-10-31T08:09:35.030Z
Answers.Text: @gc5 your awk command could have been replaced with `bedtools slop` [ref](https://bedtools.readthedocs.io/en/latest/content/tools/slop.html) which has the added advantage of making sure the intervals do not extend more than the chromosome length. It also reduces the chances of zero or one base coordinate errors (tried and test tool).
Answers.UserDisplayName: fridaymeetssunday
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-10-31T17:02:20.453Z
Answers.Text: @fridaymeetssunday amazing! I didn't know it, thanks
Answers.UserDisplayName: gc5
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-10-28T00:42:54.767Z
Answers.LastActivityDate: 2018-11-03T17:56:09.187Z
Answers.LastEditDate: 2018-11-03T17:56:09.187Z
Answers.UserDisplayName: gc5
Body: <p>I would like to filter a bed file with intervals, ie. in the format of:</p>
  
  <pre><code>chr1    13800   14301
  chr1    15500   16001
  chr1    19400   19901
  chr1    22800   23301
  </code></pre>
  
  <p>In particular, I want to filter out intervals that fall far from known genes, let's say > 1KB.</p>
  
  <p>I was thinking of querying each interval (in parallel) and check through a known set of coordinates from a gtf, e.g. from <code>Homo_sapiens.GRCh38.93.gtf</code> by filtering for gene coordinates:</p>
  
  <pre><code>1       havana  gene    168317497       168376876       .       +       .       gene_id "ENSG00000227722"; gene_version "2"; gene_name "AL022100.1"; gene_source "havana"; gene_biotype "transcribed_processed_pseudogene";
  1       mirbase gene    168375524       168375621       .       +       .       gene_id "ENSG00000207974"; gene_version "1"; gene_name "MIR557"; gene_source "mirbase"; gene_biotype "miRNA";
  1       havana  gene    168400829       168495685       .       -       .       gene_id "ENSG00000228697"; gene_version "2"; gene_name "AL023755.1"; gene_source "havana"; gene_biotype "lincRNA";
  1       havana  gene    168401483       168407274       .       -       .       gene_id "ENSG00000283255"; gene_version "1"; gene_name "AL022100.2"; gene_source "havana"; gene_biotype "lincRNA";
  1       havana  gene    168449672       168451211       .       -       .       gene_id "ENSG00000215833"; gene_version "3"; gene_name "QRSL1P1"; gene_source "havana"; gene_biotype "processed_pseudogene";
  </code></pre>
  
  <p>However, since the bedfile is quite big I am looking for the fastest solution since checking each interval may be quite expensive.</p>
  
Tags: <gene><bed><gtf><interval>
Title: Fast filtering of intervals not falling within a certain distance from known genes
ThreadText: # Fast filtering of intervals not falling within a certain distance from known genes
  
  (Asked by: username_0 on 2018-10-28)
  
  I would like to filter a bed file with intervals, ie. in the format of:
  
  ```
  chr1    13800   14301
  chr1    15500   16001
  chr1    19400   19901
  chr1    22800   23301
  
  ```
  
  In particular, I want to filter out intervals that fall far from known genes, let's say > 1KB.
  
  I was thinking of querying each interval (in parallel) and check through a known set of coordinates from a gtf, e.g. from `Homo_sapiens.GRCh38.93.gtf` by filtering for gene coordinates:
  
  ```
  1       havana  gene    168317497       168376876       .       +       .       gene_id "ENSG00000227722"; gene_version "2"; gene_name "AL022100.1"; gene_source "havana"; gene_biotype "transcribed_processed_pseudogene";
  1       mirbase gene    168375524       168375621       .       +       .       gene_id "ENSG00000207974"; gene_version "1"; gene_name "MIR557"; gene_source "mirbase"; gene_biotype "miRNA";
  1       havana  gene    168400829       168495685       .       -       .       gene_id "ENSG00000228697"; gene_version "2"; gene_name "AL023755.1"; gene_source "havana"; gene_biotype "lincRNA";
  1       havana  gene    168401483       168407274       .       -       .       gene_id "ENSG00000283255"; gene_version "1"; gene_name "AL022100.2"; gene_source "havana"; gene_biotype "lincRNA";
  1       havana  gene    168449672       168451211       .       -       .       gene_id "ENSG00000215833"; gene_version "3"; gene_name "QRSL1P1"; gene_source "havana"; gene_biotype "processed_pseudogene";
  
  ```
  
  However, since the bedfile is quite big I am looking for the fastest solution since checking each interval may be quite expensive.
  
  **Tags:** gene, bed, gtf, interval
  
  
  ---
  
  ## 3 answers
  
  ### Answer by username_1 on 2018-10-29. Score: 5
  
  Here's a way to use [BEDOPS](https://github.com/bedops/bedops), which was designed to work fast by using sorted input. Other tools now use sorting to accomplish similar performance benefits.
  
  1. Convert GTF annotations to a sorted BED file of genes:
  
     ```
     $ awk '($3=="gene")' annotations.gtf | gtf2bed - > genes.bed
  
     ```
  2. Sort your intervals, if unsorted:
  
     ```
     $ sort-bed intervals.unsorted.bed > intervals.bed
  
     ```
  3. Filter intervals against 1kb-padded genes:
  
     ```
     $ bedops --merge genes.bed | bedops --range 1000 --everything - | bedops --not-element-of 1 intervals.bed - > intervals.filtered.bed
  
     ```
  
  We first `--merge` genes to simplify dealing with overlapping genes and further speed downstream filtering.
  
  The file `intervals.filtered.bed` will contain intervals from `intervals.bed` which lie outside of a 1kb-window around merged gene space.
  
  ---
  
  To run all of this as a one-liner (assuming use of the `bash` shell):
  
  ```
  $ bedops --merge <(awk '($3=="gene")' annotations.gtf | gtf2bed -) | bedops --range 1000 --everything - | bedops --not-element-of 1 <(sort-bed intervals.unsorted.bed) - > intervals.filtered.bed
  
  ```
  
  This creates the same output result, but avoids creation of intermediate files and so should run even faster.
  
  ---
  
  **Note**: BEDOPS is agnostic about chromosome sizes — this toolkit works with any genome assembly. However, using `--range` takes care of ranges that would have a negative start position, setting their start coordinates to `0`. Bounds aren't really relevant here, in any case, because the `bedops --not-element-of` step queries input intervals, which should already be within chromosome bounds, unless they were generated incorrectly.
  
  > **Comments:**
  > **username_0**: I think there's a missing quote on the first `awk` call  
  > **username_0**: When executing 1 I get `Error: Potentially missing gene or transcript ID from GTF attributes (malformed GTF at line [1]?)`  
  > **username_1**: https://www.biostars.org/p/206342/#206350  
  
  ---
  
  ### Answer by username_2 on 2018-11-01. Score: 5
  
  In one line, using `bedtools`
  
  ```
  zcat  Homo_sapiens.GRCh38.93.gtf.gz \
   | awk '$3=="gene"' \
   | bedtools slop -b 10000 -g contigs.tsv -i - \
   | bedtools intersect -u -a intervals.bed -b - 
  
  ```
  
  This first takes the genes and filters them for the third column being gene. This is important because all `GTF` lines contain the word `gene` as `gene_id` is a manditory attribute for `GTF`.
  
  Next it passes it to slop, which adds 1 kb either side (ensuring not to go below 0 or off the end of the chr: this isn't that important because we are never outputting these regions anyway, but `bedtools` requires it). We tell `bedtools` to work on the standard in by using `-i -`.
  
  Finally we intersect these gene regions with the `intervals.bed`. `-u` says output the line from `-a` if it overlaps with at least one entry in `-b` and passing `-b -` means read file b from the standard in.
  
  This will only work if your `GTF` contains `gene` lines, which the GENCODE/Ensembl ones do. However, it is interesting to note that these GTFs are not actaully valid GTF files, as to be valid a GTF file must have both a `gene_id` and `transcript_id` attribute, and the gene lines in the GENCODE GTFs to not have this.
  
  If your GTF does not contain gene lines, but can create them by merging all the exons in all the transcripts and filling in the gaps using `gtf2gtf` from `cgat`.
  
  ```
   cgat gtf2gtf --method=merge-transcripts --stdin Homo_sapiens.GRCh38.93.gtf.gz \
   | bedtools slop -b 10000 -g contigs.tsv -i - \
   | bedtools intersect -u -a intervals.bed -b -
  
  ```
  
  
  ---
  
  ### Answer by username_0 on 2018-10-28. Score: 3
  
  First I prepared a bed file in which the gene intervals are augmented by 1KB before and after the gene start and end coordinates. Then I intersected this bed file with my original one with the option `-wa`, therefore retaining only the intervals in my original bed file that intersect with the bedfile produced from the original gtf, filtered by gene regions. The `-u` flag is also used to print only one line per overlap between the two files, in case an interval in the first overlaps with multiple intervals in the second.
  
  ```
  # download information about chromosomes
  wget http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.chrom.sizes ./
  
  # Original GTF has a different notation than UCSC
  cat Homo_sapiens.GRCh38.93.gtf | \
    sed -e $'s/^./chr&/' -e $'s/^chrMT/chrM/' | \
    awk '$3=="gene"' | \
    bedtools slop -b 1000 -g hg38.chrom.sizes -i - | \
    sed -e $'s/\tgene_id.*gene_name/\t/' -e $'s/gene_source.*//' | \
    tr -d '"; ' | \
    awk -v OFS=$'\t' '{print $1, $4, $5, $9}' > genes_1kb.bed
  
  bedtools intersect -wa -u -a intervals.bed -b genes_1kb.bed > intervals.1kb.bed
  
  ```
  
  > **Comments:**
  > **username_2**: note that this will filter out intervals that overlap a gene, or 1kp upstream, but will leave features that are within 1kbp downstream. And it will only work if the GTF has gene lines, which many GTF files don't have.  
  > **username_0**: @username_2, why this will filter out intervals that overlap a gne or 1kb upstream?  
  > **username_0**: Oops I noticed now I forgot to add 1kb to the end. Thanks! Does it solve the problem?  
  > **username_3**: @username_0 your awk command could have been replaced with `bedtools slop` [ref](https://bedtools.readthedocs.io/en/latest/content/tools/slop.html) which has the added advantage of making sure the intervals do not extend more than the chromosome length. It also reduces the chances of zero or one base coordinate errors (tried and test tool).  
  > **username_0**: @username_3 amazing! I didn't know it, thanks  
  
  ---
  
  