Entering edit mode
2.3 years ago
anasjamshed
▴
130
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
cwd
is, then look there.... Or provide absolute path for yourmy_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: