How to identify the bp distance of blast results to closest genes in a genome
2
0
Entering edit mode
4.3 years ago
FatihSarigol ▴ 260

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!

blast genome sequence • 1.4k views
ADD COMMENT
5
Entering edit mode
4.3 years ago

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

ADD COMMENT
1
Entering edit mode

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!

ADD REPLY
0
Entering edit mode

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
ADD REPLY
1
Entering edit mode
4.3 years ago

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

ADD COMMENT

Login before adding your answer.

Traffic: 1595 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6