Entering edit mode
10 months ago
anasjamshed ▴ 120
I have an unknown gene segment in the Human_gene.txt file and I want to run blastx (translated nucleotide) using the blast module of Biopython by making the E-value threshold 0.0001 and displaying the match result of 50 residues of query and subject.
I am trying this code:
import Bio from Bio.Blast import NCBIWWW from Bio import SeqIO string = open("C://Users//Home//Desktop//tsagaye work//Human_gene.txt").read() result_handle = NCBIWWW.qblast("blastx", "nr", string) # store results in a file save_file = open("my_blast.xml", "w") save_file.write(result_handle.read()) save_file.close() result_handle.close() # get BLAST results from file result_handle = open("my_blast.xml") from Bio.Blast import NCBIXML blast_record = NCBIXML.read(result_handle) #File to save results save_file = open("blastx.fasta", "w") # inspect results for alignment in blast_record.alignments: for hsp in alignment.hsps: if hsp.expect < 0.0001: print("\n" + alignment.title) print(hsp) save_file.write('%s%s%s%s' % (alignment.title, hsp.expect, alignment.length, hsp.sbjct))
But the problem is that it's creating an empty result file(blastx. fasta). Also, I want to know how can I display 50 records?
Hi, there can be multiple problems. Start with finding out, if you have any HSPs in the "my_blast.xml"... (runnning the BLAST query online can help - you can also download the BLAST xml output there and debug the result parsing portion of the code with it..)
I want to do it through code only
Ok, then inspect the xml file...
why it does not save the file?
it should (unless you are deleting it somewhere else in the code...)
Find out what your
cwdis, then look there.... Or provide absolute path for your
my_blast.xml(same way as you do with 'human_gene' file).
ok, but how can I just show 50 records?
what do you consider as your "record", hit or HSP?
you can limit it like: