1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
|
# Copyright 2005 by Michiel de Hoon. All rights reserved.
# This code is part of the Biopython distribution and governed by its
# license. Please see the LICENSE file that should have been included
# as part of this package.
import os
from Bio.Blast import NCBIXML
all_tests = [
'xbt001.xml', # BLASTP 2.2.12, gi|49176427|ref|NP_418280.3|
'xbt002.xml', # BLASTN 2.2.12, gi|1348916|gb|G26684.1|G26684
'xbt003.xml', # BLASTX 2.2.12, gi|1347369|gb|G25137.1|G25137
'xbt004.xml', # TBLASTN 2.2.12, gi|729325|sp|P39483|DHG2_BACME
'xbt005.xml', # TBLASTX 2.2.12, gi|1348853|gb|G26621.1|G26621, BLOSUM80
]
detailed_tests = [
'xbt001.xml', # BLASTP 2.2.12, gi|49176427|ref|NP_418280.3|
'xbt002.xml', # BLASTN 2.2.12, gi|1348916|gb|G26684.1|G26684
'xbt003.xml', # BLASTX 2.2.12, gi|1347369|gb|G25137.1|G25137
]
### NCBIXML.BlastParser
print "Running tests on NCBIXML.BlastParser"
for test in all_tests:
print "*" * 50, "TESTING %s" % test
datafile = os.path.join("Blast", test)
input = open(datafile)
records = NCBIXML.parse(input)
if not test in detailed_tests:
# Pull out the record, but don't print anything
for record in records:
pass
continue
for record in records:
alignments = record.alignments
print 'Blast of', record.query
print 'Found %s alignments with a total of %s HSPs' % (len(alignments),
reduce(lambda a,b: a+b, [len(a.hsps) for a in alignments]))
for alignment in alignments:
print alignment.title[:50], alignment.length, 'bp', len(alignment.hsps), 'HSPs'
# Cookbook example
E_VALUE_THRESH = 0.04
for alignment in alignments:
for hsp in alignment.hsps:
if hsp.expect < E_VALUE_THRESH:
print '*****'
print 'sequence', alignment.title
print 'length', alignment.length
# The following is a quick and dirty hack to get the
# number of digits in the exponent to match that on record
# as the expected output.
f = str(hsp.expect)
if f.find("e-") <> -1 :
matissa, exponent = str(f).split("e-")
print 'e value %se-%02i' % (matissa, int(exponent))
else :
print 'e value', hsp.expect
print hsp.query[:75] + '...'
print hsp.match[:75] + '...'
print hsp.sbjct[:75] + '...'
|
