Back to index

python-biopython  1.60
local_blast.py
Go to the documentation of this file.
00001 #!/usr/bin/env python
00002 """Script demonstrating the ability to interact with local BLAST.
00003 
00004 The contents of this script are described more fully in the available
00005 documentation.
00006 """
00007 # standard library
00008 import os
00009 import sys
00010 
00011 # biopython
00012 from Bio.Blast import NCBIStandalone
00013 
00014 my_blast_db = os.path.join(os.getcwd(), 'at-est', 'a_cds-10-7.fasta')
00015 my_blast_file = os.path.join(os.getcwd(), 'at-est', 'test_blast',
00016                              'sorghum_est-test.fasta')
00017 my_blast_exe = os.path.join(os.getcwd(), 'blast', 'blastall')
00018 
00019 print 'Running blastall...'
00020 blast_out, error_info = NCBIStandalone.blastall(my_blast_exe, 'blastn',
00021                                                 my_blast_db, my_blast_file)
00022 
00023 
00024 b_parser = NCBIStandalone.BlastParser()
00025 
00026 b_iterator = NCBIStandalone.Iterator(blast_out, b_parser)
00027 
00028 while 1:
00029     b_record = b_iterator.next()
00030 
00031     if b_record is None:
00032         break
00033     
00034     E_VALUE_THRESH = 0.04
00035     for alignment in b_record.alignments:
00036         for hsp in alignment.hsps:
00037             if hsp.expect < E_VALUE_THRESH:
00038                 print '****Alignment****'
00039                 print 'sequence:', alignment.title
00040                 print 'length:', alignment.length
00041                 print 'e value:', hsp.expect
00042 
00043                 if len(hsp.query) > 75:
00044                     dots = '...'
00045                 else:
00046                     dots = ''
00047                 
00048                 print hsp.query[0:75] + dots
00049                 print hsp.match[0:75] + dots
00050                 print hsp.sbjct[0:75] + dots
00051 
00052