BlastX through Biopython
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?
• 188 views
Traffic: 2383 users visited in the last hour
Read more here: Source link