BlastX through Biopython

Question

[anasjamshed]

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?

[massa.kassa.sc3na]

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..)

[anasjamshed]

I want to do it through code only

[massa.kassa.sc3na]

Ok, then inspect the xml file...

[anasjamshed]

why it does not save the file?

[massa.kassa.sc3na]

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 your my_blast.xml (same way as you do with 'human_gene' file).

[anasjamshed]

ok, but how can I just show 50 records?

[massa.kassa.sc3na]

what do you consider as your "record", hit or HSP?

[anasjamshed]

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))

[massa.kassa.sc3na]

you can limit it like:

iterations = 50
i=0
for a in b:
  for c in a:
    if smt:
      i += 1
    if i >= iterations:
      break
  if i >=iterations:
    break