Ok, so I am using the NcbiblastxCommandline wrapper to execute a BLAST search on a custom database.
cline = NcbiblastxCommandline(cmd='blastx', query=temp_path, out=blast_path, subject=path, outfmt=5, max_target_seqs=1)
os.system(str(cline))
The funny thing is that it is working just fine when I run those from the command line. However, inside my script, when I execute the exact same commands, I end up getting a blank xml results file instead of the actual results. I thought that maybe the script was not waiting for the process to complete, so I tried using subprocess and subprocess.wait(). That didn't work. I also tried using something else like time.sleep(60) which is plenty time enough for it to finish. That didn't work. So I am a bit out of ideas and am very puzzled as to why it works perfectly when executing as one-liners in the command shell after >>python but won't work from a script. Thanks for any suggestions!
I'm can't be certain, but the cmd='blastx' is extra to how I am used to calling it, and you don't specify a database in your call (implied?)
This is how I run it:
def blastn(fas, db):
print("blasting against the " + db + " database")
blast_in = (fas)
blast_out = (tmpdir + "blast_files/" + sample + ".blast")
blastn_cline = NcbiblastnCommandline(query=blast_in, db=db, evalue=0.001, outfmt=6, out=blast_out, num_threads=12)
stdout, stderr = blastn_cline()
blast_err_log = open(tmpdir + "blast_err.txt", "w")
blast_stdout_log = open(tmpdir + "blast_stdout.txt", "w")
blast_err_log.write(stderr)
blast_stdout_log.write(stdout)
return blast_out
Hope this helps.
Yes, the command was just copied from your post - what I'm interested in is how this looks:
stdout, stderr = blastn_cline()
print(stdout)
print(stderr)
Does stderr say anything special?
@Mabeuf So the DB is specified with the subject=path part of the expression. You use subject when you want to use a custom database, instead of one like nr etc... There are no errors and no stdout, gives an exit status of 0. I'll try it with a function definition and see if that changes anything. As I said, it runs perfectly well from the python interpreter via command line but not when placed inside a script. It does sit in a loop in my program, so maybe that has something to do with it. I really am out of ideas.
Personally, I've never used the 'subject=' call, I use 'db=' even for custom databases, so I can't comment on that. If you look at http://biopython.org/DIST/docs/api/Bio.Blast.Applications-pysrc.html at line 651 it says that "subject" is for "Subject sequence(s) to search." rather than a path to db, although it's not in the 'NcbiblastnCommandline' function. I'm just guessing though. Maybe worth a go?
No, you still use -db with a custom database (created using the makeblastdb command). You use -subject if you don't have a database, and instead want to search against a FASTA file doing pairwise comparisons (which will give potentially misleading e-values). See http://blastedbio.blogspot.co.uk/2012/05/blast-ingoring-search-space-size-for-e.html
Oh, ok. Thanks for the tip. So I went ahead and created a BLAST DB using makeblastdb which created 3 files with the extensions .phr .pin .psq (all begin with blastdb). So then I tried to rerun the blastx command using the -db flag and pointing to the blastdb. I ended up getting an xml file with no hits in it. Blast ran smoothly i.e. it found the newly created database but just didn't yield any hits in the xml file. Upon examining the .phr .pin .psq files, they have some very strange characters in them like
UPDATE: So i fixed the BLAST output problem when running it from the command line. I was using blastx instead of blastp and this, I think, was incompatible with the database type I was constructing with makeblastdb (type 'prot'). But I still cant get it to run properly from a Python program... still get the blank XML document when I should get a decent sized amount of hits. Any suggestions?
Log in to answer this question.
Have you tried to run it like this?
What do
stdoutandstderrsay?In addition to Philipp's good question, add a
print(str(cline))to see what the command Biopython tries to run for you is exactly.If you really want to use
os.system(...)then you need to check for a non-zero return code (indicating an error).