Fast filtering of intervals not falling within a certain distance from known genes

Question

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.

Answer 1

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.

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.

Answer 2

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