<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>

CC BY-SA 4.0
2018-11-01T19:16:54.723Z
I think there's a missing quote on the first `awk` call
gc5
CC BY-SA 4.0
2018-11-02T21:40:23.507Z
When executing 1 I get `Error: Potentially missing gene or transcript ID from GTF attributes (malformed GTF at line [1]?)`
gc5
CC BY-SA 4.0
2018-11-02T23:01:28.997Z
https://www.biostars.org/p/206342/#206350
Alex Reynolds
CC BY-SA 4.0
2018-10-29T17:08:16.627Z
2018-11-02T14:42:35.780Z
2018-11-02T14:42:35.780Z
Alex Reynolds
<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>

CC BY-SA 4.0
2018-11-01T19:58:53.137Z
2018-11-04T12:51:52.710Z
2018-11-04T12:51:52.710Z
Ian Sudbery
<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>

CC BY-SA 4.0
2018-10-29T12:26:00.107Z
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.
Ian Sudbery
CC BY-SA 4.0
2018-10-30T21:18:34.110Z
@IanSudbery, why this will filter out intervals that overlap a gne or 1kb upstream?
gc5
CC BY-SA 4.0
2018-10-30T22:38:38.127Z
Oops I noticed now I forgot to add 1kb to the end. Thanks! Does it solve the problem?
gc5
CC BY-SA 4.0
2018-10-31T08:09:35.030Z
@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).
fridaymeetssunday
CC BY-SA 4.0
2018-10-31T17:02:20.453Z
@fridaymeetssunday amazing! I didn't know it, thanks
gc5
CC BY-SA 4.0
2018-10-28T00:42:54.767Z
2018-11-03T17:56:09.187Z
2018-11-03T17:56:09.187Z
gc5
<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>

<gene><bed><gtf><interval>
Fast filtering of intervals not falling within a certain distance from known genes
# 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, which was designed to work fast by using sorted input. Other tools now use sorting to accomplish similar performance benefits.

Convert GTF annotations to a sorted BED file of genes:
$ awk '($3=="gene")' annotations.gtf | gtf2bed - > genes.bed

Sort your intervals, if unsorted:
$ sort-bed intervals.unsorted.bed > intervals.bed

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**: @IanSudbery, 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**: @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).  
> **username_0**: @fridaymeetssunday amazing! I didn't know it, thanks  

---

