How to obtain 200 bp flanking sequence from both sides of BLAST search?

Question

[anikcropscience]

Hello,

I have created a local blast database and blasted some genes against the reference using

blastn -query a.fasta -db BLASTDATABASE -max_target_seqs 5 -outfmt 6 -out output

I obtain a ~98% similarity score for the best hit. Now, I want to extract the fasta sequence of that BLAST search together with the 200 bp flanking sequence from both start and stop sides for designing a primer.

Can you please give me a suggestion on how to do this for a single BLAST search?

Thank you.

[manaswwm]

One way to do this, given that you have start and stop coordinates of the best hit, is by using samtools.

samtools faidx a.fasta Chr:start-stop > seq.fa

here a.fasta - reference file using in blastn; Chr - chromosome of the hit sequence; for start and stop you can add 200bp to the coordinates of the best hit and extract it from the genome

[lieven.sterck]

here you're extracting context for the query, not the hit.

[anikcropscience]

Ok. Thank you. I will try this and see if it works.

[lieven.sterck]

you want to extract context region for the hit or the query ?

and spoiler alert: you can't so this with a single BLAST search! (best will be two-step approach)

[anikcropscience]

So, I want to extract the context region for the hit. Because the query gene is from another species and I am trying to find the homologous sequence in my reference species.

Answer 1

[lieven.sterck]

For extracting the region + context of the hit, you're best of with using blastdbcmd. it has an option -range that you can use to extract a certain region from a hit in the blastdb .

As with the approach manaswwm proposed (for the query) you can extract the start & stop from the blast hit table and use that in the blastdbcmd (potentially recalculating them with + and - 200 bp )

[anikcropscience]

Thanks a lot for your helpful comments.