blast fasta file

fastablast

#! /usr/local/bin/python
#modified from Bienvenido Velez
import sys
import os.path
from Bio import SeqIO
from Bio.Blast import NCBIXML
from Bio.Blast import NCBIWWW

def runBlastn(fastaSequence,outputFileName):
    'Conducts a nucleotide Blast run on the given fastaSequence and save the blast run output in XML format'
    result_handle = NCBIWWW.qblast("blastn", "nr",fastaSequence)
    # save the blast results for later, in case we want to look at them
    save_file = open(outputFileName, "w")
    blast_results = result_handle.read()
    save_file.write(blast_results)
    save_file.close()
    return blast_results    

###########################
all_contigs_file = open("contigs.fa", "rU")
allcontigs = list(SeqIO.parse(all_contigs_file, "fasta"))
outfile = open(sys.argv[1]+'.blastresults','w')

for contig in allcontigs:
    print contig.name
    outfile.write(str(contig.name))
    if os.path.isfile(contig.name+'.out') == False:
        print 'already blast '+contig.name
        blast_results=runBlastn(contig.seq,contig.name+'.out')
    blast_out = open(contig.name+'.out')
    blast_record = NCBIXML.read(blast_out)
    
    for alignment in blast_record.alignments[0:2]:
        for hsp in alignment.hsps:
            if hsp.expect < 0.1:
                outfile.write(str(alignment.title))
                print alignment.title

outfile.close()