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!
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
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!
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
What about using bedtools intersect with increasing intervals (== expand the up/downstream region of the genes gradually)?
Log in to answer this question.
