This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Extract subject sequences from Blast result

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.

sequence blast homology
• 2,534 views • link
written by mnsp088 • mnsp088

1 answer

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

• 0 views • link
written by Nitin Narwade Nitin Narwade

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"
• 0 views • link
written by mnsp088 • mnsp088

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
• 0 views • link
written by Nitin Narwade Nitin Narwade

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

• 0 views • link
written by mnsp088 • mnsp088

Log in to answer this question.

More posts like this