Extract subject sequences from Blast result

Question

[mnsp088]

Hi all,

I ran TBLASTN (query=protein, subject=nucleotide) and got the following result:

qseqid       sseqid     pident      sstart  send    evalue  
MP_11765    contig_17   95.322      5460    5221    0

I noticed that sstart > ssend meaning this query is mapping on the reverse complement of the subject. In this case how can I pull out the DNA sequence from my subject nucleotide file that corresponds to this hit? Any ideas on a program that can do this?

Thank you.

Answer 1

[Nitin Narwade]

bedtools getfasta utility may help you. You need to include strand information in your blast output and use -s option of bedtools getfasta.

[mnsp088]

I tried pulling out strand info by adding ''sstrand" in my tblastn command below however this just results in a column of N/As

tblastn -query prot_seqs.txt -db my_nuc_db -outfmt "6 qseqid sseqid pident length qlen qstart qend sstart send evalue sstrand sseq" 

[Nitin Narwade]

Hmmmmmm!!

I am not sure why this is happening but you can add strand information manually for the blast output. Simple awk should work. Something like this,

awk -F'\t' '{ if($4 > $5) {print $2"\t"$4"\t"$5"\t-"} else{print $2"\t"$4"\t"$5"\t+"} }' blast.output.tsv >subject-with-strand-info.bed

[mnsp088]

I didn't think about adding them in manually but this works! Thank you!