Biopython Ncbiblastncommandline() Does Not Output Xml File.

Question

[chris.richard.rivera]

Hello everyone, I am new to 'real' bioinformatics and have been trying to call blastn with the biopython NcbiblastnCommandline() function on a local fasta database. It does not return an error and fails to output a xml file. My blast programs work and are in the path. I can blast against the database from the terminal and I am able to access the web based NCBI blast using biopython. My os is OS X 10.8.3.

Example code:

import os 
from Bio import SeqIO
from Bio.Blast.Applications import NcbiblastnCommandline as blastn

seqs = [seq_record for seq_record in SeqIO.parse('spiking.fasta','fasta')]
os.chdir('/Users/christopherrichardrivera/Desktop/blastdb')
blast = blastn(query=TheSeq, db='repSetBlast', outfmt=5, out='blastresult.xml')

There is no output to screen, printing blast returns:

blastn -out sequence1.xml -outfmt 5 -query ID: AB

Name: AB

Description: AB

Number of features: 0

Seq('CCTACGGGAGGCAGCAGTGGGGAATATTGGACAATGGGGGGAACCCTGATCCAG...AAT', SingleLetterAlphabet()) -db repSet -evalue 0.01

Thank you very much.

Answer 1

[Peter]

Unless you've omitted a line (and we can't tell as the code is incomplete), you never actually ran BLAST, all you did was setup the command ready to go. Try:

stdout, stderr = blast()

Does that give anything (eg an error saying command not found)?

Also note the query argument to blastn should be a filename (in FASTA format), and not a SeqRecord object etc.

[chris.richard.rivera]

Thanks for the help. I tried your recommendations, and it seems that using the line :

stdout, stderr = blast()

resulted in execution of the command and creation of the xml file.
When I simply call blast() it does not appear to create the xml file. Is it necessary to always have the stdout, stderr for this to work? If so why? Thanks

[Peter]

These should both run BLAST,

stdout, stderr = blast()

and

blast()

The second version just discards the output, the first version stores it in the two given variables.

Answer 2

[Clint Valentine]

You can also use this this line which will create an .xml output file in the current working directory.

blastn(query=TheSeq, db='repSetBlast', outfmt=5, out='blastresult.xml')()

instead of this

blast = blastn(query=TheSeq, db='repSetBlast', outfmt=5, out='blastresult.xml')
stdout, stderr = blast()