BlastX through Biopython
0
0
Entering edit mode
2.0 years 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?

Biopython • 1.3k views
ADD COMMENT
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

I want to do it through code only

ADD REPLY
0
Entering edit mode

Ok, then inspect the xml file...

ADD REPLY
0
Entering edit mode

why it does not save the file?

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

ok, but how can I just show 50 records?

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode
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))
ADD REPLY
0
Entering edit mode

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
ADD REPLY

Login before adding your answer.

Traffic: 1918 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6