Back to index

python-biopython  1.60
www_blast.py
Go to the documentation of this file.
00001 #!/usr/bin/env python
00002 """Example showing how to deal with internet BLAST from Biopython.
00003 
00004 This code is described in great detail in the BLAST section of the Biopython
00005 documentation.
00006 """
00007 # standard library
00008 import cStringIO
00009 
00010 # biopython
00011 from Bio.Blast import NCBIWWW
00012 from Bio import Fasta
00013 
00014 # first get the sequence we want to parse from a FASTA file
00015 file_for_blast = open('m_cold.fasta', 'r')
00016 f_iterator = Fasta.Iterator(file_for_blast)
00017 
00018 f_record = f_iterator.next()
00019 
00020 print 'Doing the BLAST and retrieving the results...'
00021 result_handle = NCBIWWW.qblast('blastn', 'nr', f_record)
00022 
00023 # save the results for later, in case we want to look at it
00024 save_file = open('m_cold_blast.out', 'w')
00025 blast_results = result_handle.read()
00026 save_file.write(blast_results)
00027 save_file.close()
00028 
00029 print 'Parsing the results and extracting info...'
00030 b_parser = NCBIWWW.BlastParser()
00031 
00032 # option 1 -- parse the string directly
00033 # b_record = b_parser.parse_str(blast_results)
00034 
00035 # option 2 -- create a handle from the string and parse it
00036 string_result_handle = cStringIO.StringIO(blast_results)
00037 b_record = b_parser.parse(string_result_handle)
00038 
00039 # now get the alignment info for all e values greater than some threshold
00040 E_VALUE_THRESH = 0.1
00041 
00042 for alignment in b_record.alignments:
00043     for hsp in alignment.hsps:
00044         if hsp.expect < E_VALUE_THRESH:
00045             print '****Alignment****'
00046             print 'sequence:', alignment.title
00047             print 'length:', alignment.length
00048             print 'e value:', hsp.expect
00049             print hsp.query[0:75] + '...'
00050             print hsp.match[0:75] + '...'
00051             print hsp.sbjct[0:75] + '...'
00052     
00053