Fast filtering of intervals not falling within a certain distance from known genes 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.
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 - &gt; genes.bed

Sort your intervals, if unsorted:

$ sort-bed intervals.unsorted.bed &gt; intervals.bed

Filter intervals against 1kb-padded genes:

$ bedops --merge genes.bed | bedops --range 1000 --everything - | bedops --not-element-of 1 intervals.bed - &gt; 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 &lt;(awk '($3=="gene")' annotations.gtf | gtf2bed -) | bedops --range 1000 --everything - | bedops --not-element-of 1 &lt;(sort-bed intervals.unsorted.bed) - &gt; 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.