How to identify the bp distance of blast results to closest genes in a genome

Question

[FatihSarigol]

Using the blastplus output blasted against a single species genome database prepared locally, for each blast hit, I want to search the GFF annotation file to ask how far away the hit is to the closest gene in this genome.

Easy to check manually for a few hits, but I need a program that can do it automatically. I am trying to find out if loci that are under selection are close to the coding regions to make inference if they may be controlling the expression of these flanking genes. Is there already a program that does this? If not, eventually I can share the code that I'll need to write myself here..

Thanks!

Answer 1

[Alex Reynolds]

BEDOPS gff2bed and closest-features work with BLAST output via -outfmt 6. You could use the --dist option to report the distance between the hit and its nearest element. Ref. https://bedops.readthedocs.io/en/latest/content/reference/set-operations/closest-features.html

[FatihSarigol]

thanks! though if you meant it works directly, probably something is odd about my blastplus version, it did require some parsing for my output format 6, but I got it to work anyway. great tool to have!

[Alex Reynolds]

This example returns sorted, six-column (stranded) BED from blastn, e.g.:

$ blastn -query query.fa -db /db/foo -outfmt 6 "sseqid sstart send sseqid evalue sstrand" | sort-bed - > blastHits.bed

To get genes out of a GFF file:

$ awk '($3 == "gene")' annotations.gff | gff2bed - > annotations.bed

The files blastHits.bed and annotations.bed can then be run through closest-features, to find the nearest gene to a BLAST hit.

$ closest-features --closest blastHits.bed annotations.bed > answer.bed

Answer 2

[lieven.sterck]

What about using bedtools intersect with increasing intervals (== expand the up/downstream region of the genes gradually)?