View Single Post
Old 12-02-2014, 10:38 AM   #1
skbrimer
Member
 
Location: OP Kansas

Join Date: Mar 2014
Posts: 53
Default Biopython blast parsing

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

Code:
  
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.

Last edited by skbrimer; 12-02-2014 at 10:50 AM.
skbrimer is offline   Reply With Quote