forked from biopython/biopython
-
Notifications
You must be signed in to change notification settings - Fork 0
/
NCBIStandalone.py
1818 lines (1554 loc) · 72 KB
/
NCBIStandalone.py
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
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# Copyright 1999-2000 by Jeffrey Chang. 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.
# Patches by Mike Poidinger to support multiple databases.
# Updated by Peter Cock in 2007 to do a better job on BLAST 2.2.15
"""Code for calling standalone BLAST and parsing plain text output (DEPRECATED).
Rather than parsing the human readable plain text BLAST output (which seems to
change with every update to BLAST), we and the NBCI recommend you parse the
XML output instead. The plain text parser in this module still works at the
time of writing, but is considered obsolete and updating it to cope with the
latest versions of BLAST is not a priority for us.
This module also provides code to work with the "legacy" standalone version of
NCBI BLAST, tools blastall, rpsblast and blastpgp via three helper functions of
the same name. These functions are very limited for dealing with the output as
files rather than handles, for which the wrappers in Bio.Blast.Applications are
preferred. Furthermore, the NCBI themselves regard these command line tools as
"legacy", and encourage using the new BLAST+ tools instead. Biopython has
wrappers for these under Bio.Blast.Applications (see the tutorial).
"""
from __future__ import print_function
from Bio import BiopythonDeprecationWarning
import warnings
warnings.warn("This module has been deprecated. Consider Bio.SearchIO for "
"parsing BLAST output instead.", BiopythonDeprecationWarning)
import os
import re
from Bio._py3k import StringIO
from Bio import File
from Bio.ParserSupport import *
from Bio.Blast import Record
from Bio.Application import _escape_filename
__docformat__ = "restructuredtext en"
class LowQualityBlastError(Exception):
"""Error caused by running a low quality sequence through BLAST.
When low quality sequences (like GenBank entries containing only
stretches of a single nucleotide) are BLASTed, they will result in
BLAST generating an error and not being able to perform the BLAST.
search. This error should be raised for the BLAST reports produced
in this case.
"""
pass
class ShortQueryBlastError(Exception):
"""Error caused by running a short query sequence through BLAST.
If the query sequence is too short, BLAST outputs warnings and errors::
Searching[blastall] WARNING: [000.000] AT1G08320: SetUpBlastSearch failed.
[blastall] ERROR: [000.000] AT1G08320: Blast:
[blastall] ERROR: [000.000] AT1G08320: Blast: Query must be at least wordsize
done
This exception is raised when that condition is detected.
"""
pass
class _Scanner(object):
"""Scan BLAST output from blastall or blastpgp.
Tested with blastall and blastpgp v2.0.10, v2.0.11
Methods:
- feed Feed data into the scanner.
"""
def feed(self, handle, consumer):
"""S.feed(handle, consumer)
Feed in a BLAST report for scanning. handle is a file-like
object that contains the BLAST report. consumer is a Consumer
object that will receive events as the report is scanned.
"""
if isinstance(handle, File.UndoHandle):
uhandle = handle
else:
uhandle = File.UndoHandle(handle)
# Try to fast-forward to the beginning of the blast report.
read_and_call_until(uhandle, consumer.noevent, contains='BLAST')
# Now scan the BLAST report.
self._scan_header(uhandle, consumer)
self._scan_rounds(uhandle, consumer)
self._scan_database_report(uhandle, consumer)
self._scan_parameters(uhandle, consumer)
def _scan_header(self, uhandle, consumer):
# BLASTP 2.0.10 [Aug-26-1999]
#
#
# Reference: Altschul, Stephen F., Thomas L. Madden, Alejandro A. Schaf
# Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. Lipman (1997),
# "Gapped BLAST and PSI-BLAST: a new generation of protein database sea
# programs", Nucleic Acids Res. 25:3389-3402.
#
# Query= test
# (140 letters)
#
# Database: sdqib40-1.35.seg.fa
# 1323 sequences; 223,339 total letters
#
# ========================================================
# This next example is from the online version of Blast,
# note there are TWO references, an RID line, and also
# the database is BEFORE the query line.
# Note there possibleuse of non-ASCII in the author names.
# ========================================================
#
# BLASTP 2.2.15 [Oct-15-2006]
# Reference: Altschul, Stephen F., Thomas L. Madden, Alejandro A. Sch??ffer,
# Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. Lipman
# (1997), "Gapped BLAST and PSI-BLAST: a new generation of
# protein database search programs", Nucleic Acids Res. 25:3389-3402.
#
# Reference: Sch??ffer, Alejandro A., L. Aravind, Thomas L. Madden, Sergei
# Shavirin, John L. Spouge, Yuri I. Wolf, Eugene V. Koonin, and
# Stephen F. Altschul (2001), "Improving the accuracy of PSI-BLAST
# protein database searches with composition-based statistics
# and other refinements", Nucleic Acids Res. 29:2994-3005.
#
# RID: 1166022616-19998-65316425856.BLASTQ1
#
#
# Database: All non-redundant GenBank CDS
# translations+PDB+SwissProt+PIR+PRF excluding environmental samples
# 4,254,166 sequences; 1,462,033,012 total letters
# Query= gi:16127998
# Length=428
#
consumer.start_header()
read_and_call(uhandle, consumer.version, contains='BLAST')
read_and_call_while(uhandle, consumer.noevent, blank=1)
# There might be a <pre> line, for qblast output.
attempt_read_and_call(uhandle, consumer.noevent, start="<pre>")
# Read the reference(s)
while attempt_read_and_call(uhandle,
consumer.reference, start='Reference'):
# References are normally multiline terminated by a blank line
# (or, based on the old code, the RID line)
while True:
line = uhandle.readline()
if is_blank_line(line):
consumer.noevent(line)
break
elif line.startswith("RID"):
break
else:
# More of the reference
consumer.reference(line)
# Deal with the optional RID: ...
read_and_call_while(uhandle, consumer.noevent, blank=1)
attempt_read_and_call(uhandle, consumer.reference, start="RID:")
read_and_call_while(uhandle, consumer.noevent, blank=1)
# blastpgp may have a reference for compositional score matrix
# adjustment (see Bug 2502):
if attempt_read_and_call(
uhandle, consumer.reference, start="Reference"):
read_and_call_until(uhandle, consumer.reference, blank=1)
read_and_call_while(uhandle, consumer.noevent, blank=1)
# blastpgp has a Reference for composition-based statistics.
if attempt_read_and_call(
uhandle, consumer.reference, start="Reference"):
read_and_call_until(uhandle, consumer.reference, blank=1)
read_and_call_while(uhandle, consumer.noevent, blank=1)
line = uhandle.peekline()
assert line.strip() != ""
assert not line.startswith("RID:")
if line.startswith("Query="):
# This is an old style query then database...
# Read the Query lines and the following blank line.
read_and_call(uhandle, consumer.query_info, start='Query=')
read_and_call_until(uhandle, consumer.query_info, blank=1)
read_and_call_while(uhandle, consumer.noevent, blank=1)
# Read the database lines and the following blank line.
read_and_call_until(uhandle, consumer.database_info, end='total letters')
read_and_call(uhandle, consumer.database_info, contains='sequences')
read_and_call_while(uhandle, consumer.noevent, blank=1)
elif line.startswith("Database:"):
# This is a new style database then query...
read_and_call_until(uhandle, consumer.database_info, end='total letters')
read_and_call(uhandle, consumer.database_info, contains='sequences')
read_and_call_while(uhandle, consumer.noevent, blank=1)
# Read the Query lines and the following blank line.
# Or, on BLAST 2.2.22+ there is no blank link - need to spot
# the "... Score E" line instead.
read_and_call(uhandle, consumer.query_info, start='Query=')
# BLAST 2.2.25+ has a blank line before Length=
read_and_call_until(uhandle, consumer.query_info, start='Length=')
while True:
line = uhandle.peekline()
if not line.strip() or "Score E" in line:
break
# It is more of the query (and its length)
read_and_call(uhandle, consumer.query_info)
read_and_call_while(uhandle, consumer.noevent, blank=1)
else:
raise ValueError("Invalid header?")
consumer.end_header()
def _scan_rounds(self, uhandle, consumer):
# Scan a bunch of rounds.
# Each round begins with either a "Searching......" line
# or a 'Score E' line followed by descriptions and alignments.
# The email server doesn't give the "Searching....." line.
# If there is no 'Searching.....' line then you'll first see a
# 'Results from round' line
while not self._eof(uhandle):
line = safe_peekline(uhandle)
if not line.startswith('Searching') and \
not line.startswith('Results from round') and \
re.search(r"Score +E", line) is None and \
'No hits found' not in line:
break
self._scan_descriptions(uhandle, consumer)
self._scan_alignments(uhandle, consumer)
def _scan_descriptions(self, uhandle, consumer):
# Searching..................................................done
# Results from round 2
#
#
# Sc
# Sequences producing significant alignments: (b
# Sequences used in model and found again:
#
# d1tde_2 3.4.1.4.4 (119-244) Thioredoxin reductase [Escherichia ...
# d1tcob_ 1.31.1.5.16 Calcineurin regulatory subunit (B-chain) [B...
# d1symb_ 1.31.1.2.2 Calcyclin (S100) [RAT (RATTUS NORVEGICUS)]
#
# Sequences not found previously or not previously below threshold:
#
# d1osa__ 1.31.1.5.11 Calmodulin [Paramecium tetraurelia]
# d1aoza3 2.5.1.3.3 (339-552) Ascorbate oxidase [zucchini (Cucurb...
#
# If PSI-BLAST, may also have:
#
# CONVERGED!
consumer.start_descriptions()
# Read 'Searching'
# This line seems to be missing in BLASTN 2.1.2 (others?)
attempt_read_and_call(uhandle, consumer.noevent, start='Searching')
# blastpgp 2.0.10 from NCBI 9/19/99 for Solaris sometimes crashes here.
# If this happens, the handle will yield no more information.
if not uhandle.peekline():
raise ValueError("Unexpected end of blast report. " +
"Looks suspiciously like a PSI-BLAST crash.")
# BLASTN 2.2.3 sometimes spews a bunch of warnings and errors here:
# Searching[blastall] WARNING: [000.000] AT1G08320: SetUpBlastSearch
# [blastall] ERROR: [000.000] AT1G08320: Blast:
# [blastall] ERROR: [000.000] AT1G08320: Blast: Query must be at leas
# done
# Reported by David Weisman.
# Check for these error lines and ignore them for now. Let
# the BlastErrorParser deal with them.
line = uhandle.peekline()
if "ERROR:" in line or line.startswith("done"):
read_and_call_while(uhandle, consumer.noevent, contains="ERROR:")
read_and_call(uhandle, consumer.noevent, start="done")
# Check to see if this is PSI-BLAST.
# If it is, the 'Searching' line will be followed by:
# (version 2.0.10)
# Searching.............................
# Results from round 2
# or (version 2.0.11)
# Searching.............................
#
#
# Results from round 2
# Skip a bunch of blank lines.
read_and_call_while(uhandle, consumer.noevent, blank=1)
# Check for the results line if it's there.
if attempt_read_and_call(uhandle, consumer.round, start='Results'):
read_and_call_while(uhandle, consumer.noevent, blank=1)
# Three things can happen here:
# 1. line contains 'Score E'
# 2. line contains "No hits found"
# 3. no descriptions
# The first one begins a bunch of descriptions. The last two
# indicates that no descriptions follow, and we should go straight
# to the alignments.
if not attempt_read_and_call(
uhandle, consumer.description_header,
has_re=re.compile(r'Score +E')):
# Either case 2 or 3. Look for "No hits found".
attempt_read_and_call(uhandle, consumer.no_hits,
contains='No hits found')
try:
read_and_call_while(uhandle, consumer.noevent, blank=1)
except ValueError as err:
if str(err) != "Unexpected end of stream.":
raise err
consumer.end_descriptions()
# Stop processing.
return
# Read the score header lines
read_and_call(uhandle, consumer.description_header,
start='Sequences producing')
# If PSI-BLAST, read the 'Sequences used in model' line.
attempt_read_and_call(uhandle, consumer.model_sequences,
start='Sequences used in model')
read_and_call_while(uhandle, consumer.noevent, blank=1)
# In BLAT, rather than a "No hits found" line, we just
# get no descriptions (and no alignments). This can be
# spotted because the next line is the database block:
if safe_peekline(uhandle).startswith(" Database:"):
consumer.end_descriptions()
# Stop processing.
return
# Read the descriptions and the following blank lines, making
# sure that there are descriptions.
if not uhandle.peekline().startswith('Sequences not found'):
read_and_call_until(uhandle, consumer.description, blank=1)
read_and_call_while(uhandle, consumer.noevent, blank=1)
# If PSI-BLAST, read the 'Sequences not found' line followed
# by more descriptions. However, I need to watch out for the
# case where there were no sequences not found previously, in
# which case there will be no more descriptions.
if attempt_read_and_call(uhandle, consumer.nonmodel_sequences,
start='Sequences not found'):
# Read the descriptions and the following blank lines.
read_and_call_while(uhandle, consumer.noevent, blank=1)
l = safe_peekline(uhandle)
# Brad -- added check for QUERY. On some PSI-BLAST outputs
# there will be a 'Sequences not found' line followed by no
# descriptions. Check for this case since the first thing you'll
# get is a blank line and then 'QUERY'
if not l.startswith('CONVERGED') and l[0] != '>' \
and not l.startswith('QUERY'):
read_and_call_until(uhandle, consumer.description, blank=1)
read_and_call_while(uhandle, consumer.noevent, blank=1)
attempt_read_and_call(uhandle, consumer.converged, start='CONVERGED')
read_and_call_while(uhandle, consumer.noevent, blank=1)
consumer.end_descriptions()
def _scan_alignments(self, uhandle, consumer):
if self._eof(uhandle):
return
# qblast inserts a helpful line here.
attempt_read_and_call(uhandle, consumer.noevent, start="ALIGNMENTS")
# First, check to see if I'm at the database report.
line = safe_peekline(uhandle)
if not line:
# EOF
return
elif line.startswith(' Database') or line.startswith("Lambda"):
return
elif line[0] == '>':
# XXX make a better check here between pairwise and masterslave
self._scan_pairwise_alignments(uhandle, consumer)
elif line.startswith('Effective'):
return
else:
# XXX put in a check to make sure I'm in a masterslave alignment
self._scan_masterslave_alignment(uhandle, consumer)
def _scan_pairwise_alignments(self, uhandle, consumer):
while not self._eof(uhandle):
line = safe_peekline(uhandle)
if line[0] != '>':
break
self._scan_one_pairwise_alignment(uhandle, consumer)
def _scan_one_pairwise_alignment(self, uhandle, consumer):
if self._eof(uhandle):
return
consumer.start_alignment()
self._scan_alignment_header(uhandle, consumer)
# Scan a bunch of score/alignment pairs.
while True:
if self._eof(uhandle):
# Shouldn't have issued that _scan_alignment_header event...
break
line = safe_peekline(uhandle)
if not line.startswith(' Score'):
break
self._scan_hsp(uhandle, consumer)
consumer.end_alignment()
def _scan_alignment_header(self, uhandle, consumer):
# >d1rip__ 2.24.7.1.1 Ribosomal S17 protein [Bacillus
# stearothermophilus]
# Length = 81
#
# Or, more recently with different white space:
#
# >gi|15799684|ref|NP_285696.1| threonine synthase ...
# gi|15829258|ref|NP_308031.1| threonine synthase
# ...
# Length=428
read_and_call(uhandle, consumer.title, start='>')
while True:
line = safe_readline(uhandle)
if line.lstrip().startswith('Length =') \
or line.lstrip().startswith('Length='):
consumer.length(line)
break
elif is_blank_line(line):
# Check to make sure I haven't missed the Length line
raise ValueError("I missed the Length in an alignment header")
consumer.title(line)
# Older versions of BLAST will have a line with some spaces.
# Version 2.0.14 (maybe 2.0.13?) and above print a true blank line.
if not attempt_read_and_call(uhandle, consumer.noevent,
start=' '):
read_and_call(uhandle, consumer.noevent, blank=1)
def _scan_hsp(self, uhandle, consumer):
consumer.start_hsp()
self._scan_hsp_header(uhandle, consumer)
self._scan_hsp_alignment(uhandle, consumer)
consumer.end_hsp()
def _scan_hsp_header(self, uhandle, consumer):
# Score = 22.7 bits (47), Expect = 2.5
# Identities = 10/36 (27%), Positives = 18/36 (49%)
# Strand = Plus / Plus
# Frame = +3
#
read_and_call(uhandle, consumer.score, start=' Score')
read_and_call(uhandle, consumer.identities, start=' Identities')
# BLASTN
attempt_read_and_call(uhandle, consumer.strand, start=' Strand')
# BLASTX, TBLASTN, TBLASTX
attempt_read_and_call(uhandle, consumer.frame, start=' Frame')
read_and_call(uhandle, consumer.noevent, blank=1)
def _scan_hsp_alignment(self, uhandle, consumer):
# Query: 11 GRGVSACA-------TCDGFFYRNQKVAVIGGGNTAVEEALYLSNIASEVHLIHRRDGF
# GRGVS+ TC Y + + V GGG+ + EE L + I R+
# Sbjct: 12 GRGVSSVVRRCIHKPTCKE--YAVKIIDVTGGGSFSAEEVQELREATLKEVDILRKVSG
#
# Query: 64 AEKILIKR 71
# I +K
# Sbjct: 70 PNIIQLKD 77
#
while True:
# Blastn adds an extra line filled with spaces before Query
attempt_read_and_call(uhandle, consumer.noevent, start=' ')
read_and_call(uhandle, consumer.query, start='Query')
read_and_call(uhandle, consumer.align, start=' ')
read_and_call(uhandle, consumer.sbjct, start='Sbjct')
try:
read_and_call_while(uhandle, consumer.noevent, blank=1)
except ValueError as err:
if str(err) != "Unexpected end of stream.":
raise err
# End of File (well, it looks like it with recent versions
# of BLAST for multiple queries after the Iterator class
# has broken up the whole file into chunks).
break
line = safe_peekline(uhandle)
# Alignment continues if I see a 'Query' or the spaces for Blastn.
if not (line.startswith('Query') or line.startswith(' ')):
break
def _scan_masterslave_alignment(self, uhandle, consumer):
consumer.start_alignment()
while True:
line = safe_readline(uhandle)
# Check to see whether I'm finished reading the alignment.
# This is indicated by 1) database section, 2) next psi-blast
# round, which can also be a 'Results from round' if no
# searching line is present
# patch by chapmanb
if line.startswith('Searching') or \
line.startswith('Results from round'):
uhandle.saveline(line)
break
elif line.startswith(' Database'):
uhandle.saveline(line)
break
elif is_blank_line(line):
consumer.noevent(line)
else:
consumer.multalign(line)
read_and_call_while(uhandle, consumer.noevent, blank=1)
consumer.end_alignment()
def _eof(self, uhandle):
try:
line = safe_peekline(uhandle)
except ValueError as err:
if str(err) != "Unexpected end of stream.":
raise err
line = ""
return not line
def _scan_database_report(self, uhandle, consumer):
# Database: sdqib40-1.35.seg.fa
# Posted date: Nov 1, 1999 4:25 PM
# Number of letters in database: 223,339
# Number of sequences in database: 1323
#
# Lambda K H
# 0.322 0.133 0.369
#
# Gapped
# Lambda K H
# 0.270 0.0470 0.230
#
##########################################
# Or, more recently Blast 2.2.15 gives less blank lines
##########################################
# Database: All non-redundant GenBank CDS translations+PDB+SwissProt+PIR+PRF excluding
# environmental samples
# Posted date: Dec 12, 2006 5:51 PM
# Number of letters in database: 667,088,753
# Number of sequences in database: 2,094,974
# Lambda K H
# 0.319 0.136 0.395
# Gapped
# Lambda K H
# 0.267 0.0410 0.140
if self._eof(uhandle):
return
consumer.start_database_report()
# Subset of the database(s) listed below
# Number of letters searched: 562,618,960
# Number of sequences searched: 228,924
if attempt_read_and_call(uhandle, consumer.noevent, start=" Subset"):
read_and_call(uhandle, consumer.noevent, contains="letters")
read_and_call(uhandle, consumer.noevent, contains="sequences")
read_and_call(uhandle, consumer.noevent, start=" ")
# Sameet Mehta reported seeing output from BLASTN 2.2.9 that
# was missing the "Database" stanza completely.
while attempt_read_and_call(uhandle, consumer.database,
start=' Database'):
# BLAT output ends abruptly here, without any of the other
# information. Check to see if this is the case. If so,
# then end the database report here gracefully.
if not uhandle.peekline().strip() \
or uhandle.peekline().startswith("BLAST"):
consumer.end_database_report()
return
# Database can span multiple lines.
read_and_call_until(uhandle, consumer.database, start=' Posted')
read_and_call(uhandle, consumer.posted_date, start=' Posted')
read_and_call(uhandle, consumer.num_letters_in_database,
start=' Number of letters')
read_and_call(uhandle, consumer.num_sequences_in_database,
start=' Number of sequences')
# There may not be a line starting with spaces...
attempt_read_and_call(uhandle, consumer.noevent, start=' ')
line = safe_readline(uhandle)
uhandle.saveline(line)
if 'Lambda' in line:
break
try:
read_and_call(uhandle, consumer.noevent, start='Lambda')
read_and_call(uhandle, consumer.ka_params)
except:
pass
# This blank line is optional:
attempt_read_and_call(uhandle, consumer.noevent, blank=1)
# not BLASTP
attempt_read_and_call(uhandle, consumer.gapped, start='Gapped')
# not TBLASTX
if attempt_read_and_call(uhandle, consumer.noevent, start='Lambda'):
read_and_call(uhandle, consumer.ka_params_gap)
# Blast 2.2.4 can sometimes skip the whole parameter section.
# Thus, I need to be careful not to read past the end of the
# file.
try:
read_and_call_while(uhandle, consumer.noevent, blank=1)
except ValueError as x:
if str(x) != "Unexpected end of stream.":
raise
consumer.end_database_report()
def _scan_parameters(self, uhandle, consumer):
# Matrix: BLOSUM62
# Gap Penalties: Existence: 11, Extension: 1
# Number of Hits to DB: 50604
# Number of Sequences: 1323
# Number of extensions: 1526
# Number of successful extensions: 6
# Number of sequences better than 10.0: 5
# Number of HSP's better than 10.0 without gapping: 5
# Number of HSP's successfully gapped in prelim test: 0
# Number of HSP's that attempted gapping in prelim test: 1
# Number of HSP's gapped (non-prelim): 5
# length of query: 140
# length of database: 223,339
# effective HSP length: 39
# effective length of query: 101
# effective length of database: 171,742
# effective search space: 17345942
# effective search space used: 17345942
# T: 11
# A: 40
# X1: 16 ( 7.4 bits)
# X2: 38 (14.8 bits)
# X3: 64 (24.9 bits)
# S1: 41 (21.9 bits)
# S2: 42 (20.8 bits)
##########################################
# Or, more recently Blast(x) 2.2.15 gives
##########################################
# Matrix: BLOSUM62
# Gap Penalties: Existence: 11, Extension: 1
# Number of Sequences: 4535438
# Number of Hits to DB: 2,588,844,100
# Number of extensions: 60427286
# Number of successful extensions: 126433
# Number of sequences better than 2.0: 30
# Number of HSP's gapped: 126387
# Number of HSP's successfully gapped: 35
# Length of query: 291
# Length of database: 1,573,298,872
# Length adjustment: 130
# Effective length of query: 161
# Effective length of database: 983,691,932
# Effective search space: 158374401052
# Effective search space used: 158374401052
# Neighboring words threshold: 12
# Window for multiple hits: 40
# X1: 16 ( 7.3 bits)
# X2: 38 (14.6 bits)
# X3: 64 (24.7 bits)
# S1: 41 (21.7 bits)
# S2: 32 (16.9 bits)
# Blast 2.2.4 can sometimes skip the whole parameter section.
# BLAT also skips the whole parameter section.
# Thus, check to make sure that the parameter section really
# exists.
if not uhandle.peekline().strip():
return
# BLASTN 2.2.9 looks like it reverses the "Number of Hits" and
# "Number of Sequences" lines.
consumer.start_parameters()
# Matrix line may be missing in BLASTN 2.2.9
attempt_read_and_call(uhandle, consumer.matrix, start='Matrix')
# not TBLASTX
attempt_read_and_call(uhandle, consumer.gap_penalties, start='Gap')
attempt_read_and_call(uhandle, consumer.num_sequences,
start='Number of Sequences')
attempt_read_and_call(uhandle, consumer.num_hits,
start='Number of Hits')
attempt_read_and_call(uhandle, consumer.num_sequences,
start='Number of Sequences')
attempt_read_and_call(uhandle, consumer.num_extends,
start='Number of extensions')
attempt_read_and_call(uhandle, consumer.num_good_extends,
start='Number of successful')
attempt_read_and_call(uhandle, consumer.num_seqs_better_e,
start='Number of sequences')
# not BLASTN, TBLASTX
if attempt_read_and_call(uhandle, consumer.hsps_no_gap,
start="Number of HSP's better"):
# BLASTN 2.2.9
if attempt_read_and_call(uhandle, consumer.noevent,
start="Number of HSP's gapped:"):
read_and_call(uhandle, consumer.noevent,
start="Number of HSP's successfully")
# This is omitted in 2.2.15
attempt_read_and_call(uhandle, consumer.noevent,
start="Number of extra gapped extensions")
else:
read_and_call(uhandle, consumer.hsps_prelim_gapped,
start="Number of HSP's successfully")
read_and_call(uhandle, consumer.hsps_prelim_gap_attempted,
start="Number of HSP's that")
read_and_call(uhandle, consumer.hsps_gapped,
start="Number of HSP's gapped")
# e.g. BLASTX 2.2.15 where the "better" line is missing
elif attempt_read_and_call(uhandle, consumer.noevent,
start="Number of HSP's gapped"):
read_and_call(uhandle, consumer.noevent,
start="Number of HSP's successfully")
# not in blastx 2.2.1
attempt_read_and_call(uhandle, consumer.query_length,
has_re=re.compile(r"[Ll]ength of query"))
# Not in BLASTX 2.2.22+
attempt_read_and_call(uhandle, consumer.database_length,
has_re=re.compile(r"[Ll]ength of \s*[Dd]atabase"))
# BLASTN 2.2.9
attempt_read_and_call(uhandle, consumer.noevent,
start="Length adjustment")
attempt_read_and_call(uhandle, consumer.effective_hsp_length,
start='effective HSP')
# Not in blastx 2.2.1
attempt_read_and_call(
uhandle, consumer.effective_query_length,
has_re=re.compile(r'[Ee]ffective length of query'))
# This is not in BLASTP 2.2.15
attempt_read_and_call(
uhandle, consumer.effective_database_length,
has_re=re.compile(r'[Ee]ffective length of \s*[Dd]atabase'))
# Not in blastx 2.2.1, added a ':' to distinguish between
# this and the 'effective search space used' line
attempt_read_and_call(
uhandle, consumer.effective_search_space,
has_re=re.compile(r'[Ee]ffective search space:'))
# Does not appear in BLASTP 2.0.5
attempt_read_and_call(
uhandle, consumer.effective_search_space_used,
has_re=re.compile(r'[Ee]ffective search space used'))
# BLASTX, TBLASTN, TBLASTX
attempt_read_and_call(uhandle, consumer.frameshift, start='frameshift')
# not in BLASTN 2.2.9
attempt_read_and_call(uhandle, consumer.threshold, start='T')
# In BLASTX 2.2.15 replaced by: "Neighboring words threshold: 12"
attempt_read_and_call(uhandle, consumer.threshold, start='Neighboring words threshold')
# not in BLASTX 2.2.15
attempt_read_and_call(uhandle, consumer.window_size, start='A')
# get this instead: "Window for multiple hits: 40"
attempt_read_and_call(uhandle, consumer.window_size, start='Window for multiple hits')
# not in BLASTX 2.2.22+
attempt_read_and_call(uhandle, consumer.dropoff_1st_pass, start='X1')
# not TBLASTN
attempt_read_and_call(uhandle, consumer.gap_x_dropoff, start='X2')
# not BLASTN, TBLASTX
attempt_read_and_call(uhandle, consumer.gap_x_dropoff_final,
start='X3')
# not TBLASTN
attempt_read_and_call(uhandle, consumer.gap_trigger, start='S1')
# not in blastx 2.2.1
# first we make sure we have additional lines to work with, if
# not then the file is done and we don't have a final S2
if not is_blank_line(uhandle.peekline(), allow_spaces=1):
read_and_call(uhandle, consumer.blast_cutoff, start='S2')
consumer.end_parameters()
class BlastParser(AbstractParser):
"""Parses BLAST data into a Record.Blast object.
"""
def __init__(self):
"""__init__(self)"""
self._scanner = _Scanner()
self._consumer = _BlastConsumer()
def parse(self, handle):
"""parse(self, handle)"""
self._scanner.feed(handle, self._consumer)
return self._consumer.data
class PSIBlastParser(AbstractParser):
"""Parses BLAST data into a Record.PSIBlast object.
"""
def __init__(self):
"""__init__(self)"""
self._scanner = _Scanner()
self._consumer = _PSIBlastConsumer()
def parse(self, handle):
"""parse(self, handle)"""
self._scanner.feed(handle, self._consumer)
return self._consumer.data
class _HeaderConsumer(object):
def start_header(self):
self._header = Record.Header()
def version(self, line):
c = line.split()
self._header.application = c[0]
self._header.version = c[1]
if len(c) > 2:
# The date is missing in the new C++ output from blastx 2.2.22+
# Just get "BLASTX 2.2.22+\n" and that's all.
self._header.date = c[2][1:-1]
def reference(self, line):
if line.startswith('Reference: '):
self._header.reference = line[11:]
else:
self._header.reference = self._header.reference + line
def query_info(self, line):
if line.startswith('Query= '):
self._header.query = line[7:].lstrip()
elif line.startswith('Length='):
# New style way to give the query length in BLAST 2.2.22+ (the C++ code)
self._header.query_letters = _safe_int(line[7:].strip())
elif not line.startswith(' '): # continuation of query_info
self._header.query = "%s%s" % (self._header.query, line)
else:
# Hope it is the old style way to give the query length:
letters, = _re_search(
r"([0-9,]+) letters", line,
"I could not find the number of letters in line\n%s" % line)
self._header.query_letters = _safe_int(letters)
def database_info(self, line):
line = line.rstrip()
if line.startswith('Database: '):
self._header.database = line[10:]
elif not line.endswith('total letters'):
if self._header.database:
# Need to include a space when merging multi line datase descr
self._header.database = self._header.database + " " + line.strip()
else:
self._header.database = line.strip()
else:
sequences, letters = _re_search(
r"([0-9,]+) sequences; ([0-9,-]+) total letters", line,
"I could not find the sequences and letters in line\n%s" % line)
self._header.database_sequences = _safe_int(sequences)
self._header.database_letters = _safe_int(letters)
def end_header(self):
# Get rid of the trailing newlines
self._header.reference = self._header.reference.rstrip()
self._header.query = self._header.query.rstrip()
class _DescriptionConsumer(object):
def start_descriptions(self):
self._descriptions = []
self._model_sequences = []
self._nonmodel_sequences = []
self._converged = 0
self._type = None
self._roundnum = None
self.__has_n = 0 # Does the description line contain an N value?
def description_header(self, line):
if line.startswith('Sequences producing'):
cols = line.split()
if cols[-1] == 'N':
self.__has_n = 1
def description(self, line):
dh = self._parse(line)
if self._type == 'model':
self._model_sequences.append(dh)
elif self._type == 'nonmodel':
self._nonmodel_sequences.append(dh)
else:
self._descriptions.append(dh)
def model_sequences(self, line):
self._type = 'model'
def nonmodel_sequences(self, line):
self._type = 'nonmodel'
def converged(self, line):
self._converged = 1
def no_hits(self, line):
pass
def round(self, line):
if not line.startswith('Results from round'):
raise ValueError("I didn't understand the round line\n%s" % line)
self._roundnum = _safe_int(line[18:].strip())
def end_descriptions(self):
pass
def _parse(self, description_line):
line = description_line # for convenience
dh = Record.Description()
# I need to separate the score and p-value from the title.
# sp|P21297|FLBT_CAUCR FLBT PROTEIN [snip] 284 7e-77
# sp|P21297|FLBT_CAUCR FLBT PROTEIN [snip] 284 7e-77 1
# special cases to handle:
# - title must be preserved exactly (including whitespaces)
# - score could be equal to e-value (not likely, but what if??)
# - sometimes there's an "N" score of '1'.
cols = line.split()
if len(cols) < 3:
raise ValueError(
"Line does not appear to contain description:\n%s" % line)
if self.__has_n:
i = line.rfind(cols[-1]) # find start of N
i = line.rfind(cols[-2], 0, i) # find start of p-value
i = line.rfind(cols[-3], 0, i) # find start of score
else:
i = line.rfind(cols[-1]) # find start of p-value
i = line.rfind(cols[-2], 0, i) # find start of score
if self.__has_n:
dh.title, dh.score, dh.e, dh.num_alignments = \
line[:i].rstrip(), cols[-3], cols[-2], cols[-1]
else:
dh.title, dh.score, dh.e, dh.num_alignments = \
line[:i].rstrip(), cols[-2], cols[-1], 1
dh.num_alignments = _safe_int(dh.num_alignments)
dh.score = _safe_int(dh.score)
dh.e = _safe_float(dh.e)
return dh
class _AlignmentConsumer(object):
# This is a little bit tricky. An alignment can either be a
# pairwise alignment or a multiple alignment. Since it's difficult
# to know a-priori which one the blast record will contain, I'm going
# to make one class that can parse both of them.
def start_alignment(self):
self._alignment = Record.Alignment()
self._multiple_alignment = Record.MultipleAlignment()
def title(self, line):
if self._alignment.title:
self._alignment.title += " "
self._alignment.title += line.strip()
def length(self, line):
# e.g. "Length = 81" or more recently, "Length=428"
parts = line.replace(" ", "").split("=")
assert len(parts) == 2, "Unrecognised format length line"
self._alignment.length = parts[1]
self._alignment.length = _safe_int(self._alignment.length)
def multalign(self, line):
# Standalone version uses 'QUERY', while WWW version uses blast_tmp.
if line.startswith('QUERY') or line.startswith('blast_tmp'):
# If this is the first line of the multiple alignment,
# then I need to figure out how the line is formatted.
# Format of line is:
# QUERY 1 acttg...gccagaggtggtttattcagtctccataagagaggggacaaacg 60
try:
name, start, seq, end = line.split()
except ValueError:
raise ValueError("I do not understand the line\n%s" % line)
self._start_index = line.index(start, len(name))
self._seq_index = line.index(seq,