Biopython blast parsing

[skbrimer]

Hello all,

I am currently trying to automate some of the work that I do for some in-house dsRNA viral sequencing. I'm pretty new to this but instead of re-inventing the wheel I'm using Biopython and so far it has work well following the tutorial. However my ultimate goal with this to be able to take my list of de novo contigs, BLAST them, pull the info I want from the blast into another flat file (csv) and be sort them into like groups.

This is what I'm doing manually and it takes forever! I'm not even sure its the right way to go but it is working so I shouldn't complain. It would be a huge help to be able to put this list together automatically so I can spend more time in closing the sequences instead of the bulk of the time just trying to figure out what they are!

So far my script is working (just on one of the contigs), I'm more confused on how to use the NCBIXML.parse() command

  
import sys
from Bio import SeqIO
from Bio.Blast import NCBIWWW
from Bio.Blast import NCBIXML

#this allows for any fasta file to be parsed using Biopython
file_name =sys.argv[1]

#reads in the list of cotigs
handle = open(file_name, "rU")
records = list(SeqIO.parse(handle,"fasta"))
handle.close()

#will blast the contigs
results_handle = NCBIWWW.qblast("blastn", "nr", records[0].seq, hitlist_size = 1)
#save_file = open("my_blast.xml","w")
save_file = open("my_blast.text","w")
save_file.write(results_handle.read())
save_file.close()
results_handle.close()

#this section will be for the parsing and exporting of the 
#data I'm looking for into a csv for easy sorting and visual inspection

I would like the csv file to contain the Seqrecord.id, Hit_id, Hit_def, Hsp_query-from, HSP_query-to, Hsp_Hit-from, Hsp-Hit-to, Hsp_score

Is this possible to do in a simple script? The only thing I want to keep from this is the csv at the end so I don't think I need to setup an actual db, am I correct?

EDIT:: corrected the code to show that I'm trying to get this information from the xml file not the text file.

[maubp]

Given you want to make a CSV file, why mess about with XML? One of the best things in BLAST+ is the pick-your-own columns in tabular or CSV output.

[skbrimer]

Thanks for the advice Peter,

To be able to do this would I find the documentation in Biopython or in NCBI? I ask because the output says I'm using the the current BLAST (2.2.30+) however in the Biopython documentation I have not found the where to do that, I'm sorry if this is a easy question but I'm very new to both programming and genomics and any guidance you could provide would be greatly appreciated.

Sean

[maubp]

At the command line,

$ blastn -help
...

You are looking for the -outfmt option, which can be simply a number, say 6 for tabular, or a quoted string like "6 std stitlesall" where you list the fields you want (std is shorthand for the default 12 columns).

You can do this via the BLAST+ wrapper object in Biopython too, but be careful about the quotes - see https://www.biostars.org/p/95911/

[skbrimer]

Thank you Peter!

I will work on getting this implimented. Since I do not have the command line BLAST tool on my machine... well right now at least, that is easy to fix. I will install it and see if I can get the output I'm looking for.

Thank you again. I will let you know how this works.

Sean

[skbrimer]

Hi Peter,

I just waned to say that I feel this will work. I am having issues with the implementation of using both the BLAST tool and Biopython, however your suggestion to use the BLAST+ tool is the correct one.

Thank you again,
Sean