forked from biopython/biopython
-
Notifications
You must be signed in to change notification settings - Fork 0
/
NCBIXML.py
675 lines (536 loc) · 23.9 KB
/
NCBIXML.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
# Copyright 2000 by Bertrand Frottier . All rights reserved.
# Revisions 2005-2006 copyright Michiel de Hoon
# Revisions 2006-2009 copyright Peter Cock
# 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.
"""Code to work with the BLAST XML output.
The BLAST XML DTD file is on the NCBI FTP site at:
ftp://ftp.ncbi.nlm.nih.gov/blast/documents/xml/NCBI_BlastOutput.dtd
"""
from __future__ import print_function
from Bio.Blast import Record
import xml.sax
from xml.sax.handler import ContentHandler
from functools import reduce
__docformat__ = "restructuredtext en"
class _XMLparser(ContentHandler):
"""Generic SAX Parser (PRIVATE).
Just a very basic SAX parser.
Redefine the methods startElement, characters and endElement.
"""
def __init__(self, debug=0):
"""Constructor.
debug - integer, amount of debug information to print
"""
self._tag = []
self._value = ''
self._debug = debug
self._debug_ignore_list = []
def _secure_name(self, name):
"""Removes 'dangerous' from tag names.
name -- name to be 'secured'
"""
# Replace '-' with '_' in XML tag names
return name.replace('-', '_')
def startElement(self, name, attr):
"""Found XML start tag.
No real need of attr, BLAST DTD doesn't use them
name -- name of the tag
attr -- tag attributes
"""
self._tag.append(name)
# Try to call a method (defined in subclasses)
method = self._secure_name('_start_' + name)
# Note could use try / except AttributeError
# BUT I found often triggered by nested errors...
if hasattr(self, method):
eval("self.%s()" % method)
if self._debug > 4:
print("NCBIXML: Parsed: " + method)
elif self._debug > 3:
# Doesn't exist (yet) and may want to warn about it
if method not in self._debug_ignore_list:
print("NCBIXML: Ignored: " + method)
self._debug_ignore_list.append(method)
# We don't care about white space in parent tags like Hsp,
# but that white space doesn't belong to child tags like Hsp_midline
if self._value.strip():
raise ValueError("What should we do with %s before the %s tag?"
% (repr(self._value), name))
self._value = ""
def characters(self, ch):
"""Found some text.
ch -- characters read
"""
self._value += ch # You don't ever get the whole string
def endElement(self, name):
"""Found XML end tag.
name -- tag name
"""
# DON'T strip any white space, we may need it e.g. the hsp-midline
# Try to call a method (defined in subclasses)
method = self._secure_name('_end_' + name)
# Note could use try / except AttributeError
# BUT I found often triggered by nested errors...
if hasattr(self, method):
eval("self.%s()" % method)
if self._debug > 2:
print("NCBIXML: Parsed: %s %s" % (method, self._value))
elif self._debug > 1:
# Doesn't exist (yet) and may want to warn about it
if method not in self._debug_ignore_list:
print("NCBIXML: Ignored: %s %s" % (method, self._value))
self._debug_ignore_list.append(method)
# Reset character buffer
self._value = ''
class BlastParser(_XMLparser):
"""Parse XML BLAST data into a Record.Blast object.
Parses XML output from BLAST (direct use discouraged).
This (now) returns a list of Blast records.
Historically it returned a single Blast record.
You are expected to use this via the parse or read functions.
All XML 'action' methods are private methods and may be:
- ``_start_TAG`` called when the start tag is found
- ``_end_TAG`` called when the end tag is found
"""
def __init__(self, debug=0):
"""Constructor.
debug - integer, amount of debug information to print
"""
# Calling superclass method
_XMLparser.__init__(self, debug)
self._parser = xml.sax.make_parser()
self._parser.setContentHandler(self)
# To avoid ValueError: unknown url type: NCBI_BlastOutput.dtd
self._parser.setFeature(xml.sax.handler.feature_validation, 0)
self._parser.setFeature(xml.sax.handler.feature_namespaces, 0)
self._parser.setFeature(xml.sax.handler.feature_external_pes, 0)
self._parser.setFeature(xml.sax.handler.feature_external_ges, 0)
self.reset()
def reset(self):
"""Reset all the data allowing reuse of the BlastParser() object."""
self._records = []
self._header = Record.Header()
self._parameters = Record.Parameters()
self._parameters.filter = None # Maybe I should update the class?
def _start_Iteration(self):
self._blast = Record.Blast()
pass
def _end_Iteration(self):
# We stored a lot of generic "top level" information
# in self._header (an object of type Record.Header)
self._blast.reference = self._header.reference
self._blast.date = self._header.date
self._blast.version = self._header.version
self._blast.database = self._header.database
self._blast.application = self._header.application
# These are required for "old" pre 2.2.14 files
# where only <BlastOutput_query-ID>, <BlastOutput_query-def>
# and <BlastOutput_query-len> were used. Now they
# are suplemented/replaced by <Iteration_query-ID>,
# <Iteration_query-def> and <Iteration_query-len>
if not hasattr(self._blast, "query") \
or not self._blast.query:
self._blast.query = self._header.query
if not hasattr(self._blast, "query_id") \
or not self._blast.query_id:
self._blast.query_id = self._header.query_id
if not hasattr(self._blast, "query_letters") \
or not self._blast.query_letters:
self._blast.query_letters = self._header.query_letters
# Hack to record the query length as both the query_letters and
# query_length properties (as in the plain text parser, see
# Bug 2176 comment 12):
self._blast.query_length = self._blast.query_letters
# Perhaps in the long term we should deprecate one, but I would
# prefer to drop query_letters - so we need a transition period
# with both.
# Hack to record the claimed database size as database_length
# (as well as in num_letters_in_database, see Bug 2176 comment 13):
self._blast.database_length = self._blast.num_letters_in_database
# TODO? Deprecate database_letters next?
# Hack to record the claimed database sequence count as database_sequences
self._blast.database_sequences = self._blast.num_sequences_in_database
# Apply the "top level" parameter information
self._blast.matrix = self._parameters.matrix
self._blast.num_seqs_better_e = self._parameters.num_seqs_better_e
self._blast.gap_penalties = self._parameters.gap_penalties
self._blast.filter = self._parameters.filter
self._blast.expect = self._parameters.expect
self._blast.sc_match = self._parameters.sc_match
self._blast.sc_mismatch = self._parameters.sc_mismatch
# Add to the list
self._records.append(self._blast)
# Clear the object (a new empty one is create in _start_Iteration)
self._blast = None
if self._debug:
print("NCBIXML: Added Blast record to results")
# Header
def _end_BlastOutput_program(self):
"""BLAST program, e.g., blastp, blastn, etc.
Save this to put on each blast record object
"""
self._header.application = self._value.upper()
def _end_BlastOutput_version(self):
"""version number and date of the BLAST engine.
e.g. "BLASTX 2.2.12 [Aug-07-2005]" but there can also be
variants like "BLASTP 2.2.18+" without the date.
Save this to put on each blast record object
"""
parts = self._value.split()
# TODO - Check the first word starts with BLAST?
# The version is the second word (field one)
self._header.version = parts[1]
# Check there is a third word (the date)
if len(parts) >= 3:
if parts[2][0] == "[" and parts[2][-1] == "]":
self._header.date = parts[2][1:-1]
else:
# Assume this is still a date, but without the
# square brackets
self._header.date = parts[2]
def _end_BlastOutput_reference(self):
"""a reference to the article describing the algorithm (PRIVATE).
Save this to put on each blast record object
"""
self._header.reference = self._value
def _end_BlastOutput_db(self):
"""the database(s) searched (PRIVATE).
Save this to put on each blast record object
"""
self._header.database = self._value
def _end_BlastOutput_query_ID(self):
"""the identifier of the query (PRIVATE).
Important in old pre 2.2.14 BLAST, for recent versions
<Iteration_query-ID> is enough
"""
self._header.query_id = self._value
def _end_BlastOutput_query_def(self):
"""the definition line of the query (PRIVATE).
Important in old pre 2.2.14 BLAST, for recent versions
<Iteration_query-def> is enough
"""
self._header.query = self._value
def _end_BlastOutput_query_len(self):
"""the length of the query (PRIVATE).
Important in old pre 2.2.14 BLAST, for recent versions
<Iteration_query-len> is enough
"""
self._header.query_letters = int(self._value)
def _end_Iteration_query_ID(self):
"""the identifier of the query (PRIVATE)."""
self._blast.query_id = self._value
def _end_Iteration_query_def(self):
"""the definition line of the query (PRIVATE)."""
self._blast.query = self._value
def _end_Iteration_query_len(self):
"""the length of the query (PRIVATE)."""
self._blast.query_letters = int(self._value)
# def _end_BlastOutput_query_seq(self):
# """the query sequence (PRIVATE)."""
# pass # XXX Missing in Record.Blast ?
# def _end_BlastOutput_iter_num(self):
# """the psi-blast iteration number (PRIVATE)."""
# pass # XXX TODO PSI
def _end_BlastOutput_hits(self):
"""hits to the database sequences, one for every sequence (PRIVATE)."""
self._blast.num_hits = int(self._value)
# def _end_BlastOutput_message(self):
# """error messages (PRIVATE)."""
# pass # XXX What to do ?
# Parameters
def _end_Parameters_matrix(self):
"""matrix used (-M on legacy BLAST) (PRIVATE)."""
self._parameters.matrix = self._value
def _end_Parameters_expect(self):
"""expect values cutoff (PRIVATE)."""
# NOTE: In old text output there was a line:
# Number of sequences better than 1.0e-004: 1
# As far as I can see, parameters.num_seqs_better_e
# would take the value of 1, and the expectation
# value was not recorded.
#
# Anyway we should NOT record this against num_seqs_better_e
self._parameters.expect = self._value
# def _end_Parameters_include(self):
# """inclusion threshold for a psi-blast iteration (-h) (PRIVATE)."""
# pass # XXX TODO PSI
def _end_Parameters_sc_match(self):
"""match score for nucleotide-nucleotide comparaison (-r) (PRIVATE)."""
self._parameters.sc_match = int(self._value)
def _end_Parameters_sc_mismatch(self):
"""mismatch penalty for nucleotide-nucleotide comparaison (-r) (PRIVATE)."""
self._parameters.sc_mismatch = int(self._value)
def _end_Parameters_gap_open(self):
"""gap existence cost (-G) (PRIVATE)."""
self._parameters.gap_penalties = int(self._value)
def _end_Parameters_gap_extend(self):
"""gap extension cose (-E) (PRIVATE)."""
self._parameters.gap_penalties = (self._parameters.gap_penalties,
int(self._value))
def _end_Parameters_filter(self):
"""filtering options (-F) (PRIVATE)."""
self._parameters.filter = self._value
# def _end_Parameters_pattern(self):
# """pattern used for phi-blast search
# """
# pass # XXX TODO PSI
# def _end_Parameters_entrez_query(self):
# """entrez query used to limit search
# """
# pass # XXX TODO PSI
# Hits
def _start_Hit(self):
self._blast.alignments.append(Record.Alignment())
self._blast.descriptions.append(Record.Description())
self._blast.multiple_alignment = []
self._hit = self._blast.alignments[-1]
self._descr = self._blast.descriptions[-1]
self._descr.num_alignments = 0
def _end_Hit(self):
# Cleanup
self._blast.multiple_alignment = None
self._hit = None
self._descr = None
def _end_Hit_id(self):
"""identifier of the database sequence (PRIVATE)."""
self._hit.hit_id = self._value
self._hit.title = self._value + ' '
def _end_Hit_def(self):
"""definition line of the database sequence (PRIVATE)."""
self._hit.hit_def = self._value
self._hit.title += self._value
self._descr.title = self._hit.title
def _end_Hit_accession(self):
"""accession of the database sequence (PRIVATE)."""
self._hit.accession = self._value
self._descr.accession = self._value
def _end_Hit_len(self):
self._hit.length = int(self._value)
# HSPs
def _start_Hsp(self):
# Note that self._start_Hit() should have been called
# to setup things like self._blast.multiple_alignment
self._hit.hsps.append(Record.HSP())
self._hsp = self._hit.hsps[-1]
self._descr.num_alignments += 1
self._blast.multiple_alignment.append(Record.MultipleAlignment())
self._mult_al = self._blast.multiple_alignment[-1]
# Hsp_num is useless
def _end_Hsp_score(self):
"""raw score of HSP (PRIVATE)."""
self._hsp.score = float(self._value)
if self._descr.score is None:
self._descr.score = float(self._value)
def _end_Hsp_bit_score(self):
"""bit score of HSP (PRIVATE)."""
self._hsp.bits = float(self._value)
if self._descr.bits is None:
self._descr.bits = float(self._value)
def _end_Hsp_evalue(self):
"""expect value of the HSP (PRIVATE)."""
self._hsp.expect = float(self._value)
if self._descr.e is None:
self._descr.e = float(self._value)
def _end_Hsp_query_from(self):
"""offset of query at the start of the alignment (one-offset) (PRIVATE)."""
self._hsp.query_start = int(self._value)
def _end_Hsp_query_to(self):
"""offset of query at the end of the alignment (one-offset) (PRIVATE)."""
self._hsp.query_end = int(self._value)
def _end_Hsp_hit_from(self):
"""offset of the database at the start of the alignment (one-offset) (PRIVATE)."""
self._hsp.sbjct_start = int(self._value)
def _end_Hsp_hit_to(self):
"""offset of the database at the end of the alignment (one-offset) (PRIVATE)."""
self._hsp.sbjct_end = int(self._value)
# def _end_Hsp_pattern_from(self):
# """start of phi-blast pattern on the query (one-offset) (PRIVATE)."""
# pass # XXX TODO PSI
# def _end_Hsp_pattern_to(self):
# """end of phi-blast pattern on the query (one-offset) (PRIVATE)."""
# pass # XXX TODO PSI
def _end_Hsp_query_frame(self):
"""frame of the query if applicable (PRIVATE)."""
self._hsp.frame = (int(self._value),)
def _end_Hsp_hit_frame(self):
"""frame of the database sequence if applicable (PRIVATE)."""
self._hsp.frame += (int(self._value),)
def _end_Hsp_identity(self):
"""number of identities in the alignment (PRIVATE)."""
self._hsp.identities = int(self._value)
def _end_Hsp_positive(self):
"""number of positive (conservative) substitutions in the alignment (PRIVATE)."""
self._hsp.positives = int(self._value)
def _end_Hsp_gaps(self):
"""number of gaps in the alignment (PRIVATE)."""
self._hsp.gaps = int(self._value)
def _end_Hsp_align_len(self):
"""length of the alignment (PRIVATE)."""
self._hsp.align_length = int(self._value)
# def _en_Hsp_density(self):
# """score density (PRIVATE)."""
# pass # XXX ???
def _end_Hsp_qseq(self):
"""alignment string for the query (PRIVATE)."""
self._hsp.query = self._value
def _end_Hsp_hseq(self):
"""alignment string for the database (PRIVATE)."""
self._hsp.sbjct = self._value
def _end_Hsp_midline(self):
"""Formatting middle line as normally seen in BLAST report (PRIVATE)."""
self._hsp.match = self._value # do NOT strip spaces!
assert len(self._hsp.match) == len(self._hsp.query)
assert len(self._hsp.match) == len(self._hsp.sbjct)
# Statistics
def _end_Statistics_db_num(self):
"""number of sequences in the database (PRIVATE)."""
self._blast.num_sequences_in_database = int(self._value)
def _end_Statistics_db_len(self):
"""number of letters in the database (PRIVATE)."""
self._blast.num_letters_in_database = int(self._value)
def _end_Statistics_hsp_len(self):
"""the effective HSP length (PRIVATE)."""
self._blast.effective_hsp_length = int(self._value)
def _end_Statistics_eff_space(self):
"""the effective search space (PRIVATE)."""
self._blast.effective_search_space = float(self._value)
def _end_Statistics_kappa(self):
"""Karlin-Altschul parameter K (PRIVATE)."""
self._blast.ka_params = float(self._value)
def _end_Statistics_lambda(self):
"""Karlin-Altschul parameter Lambda (PRIVATE)."""
self._blast.ka_params = (float(self._value), self._blast.ka_params)
def _end_Statistics_entropy(self):
"""Karlin-Altschul parameter H (PRIVATE)."""
self._blast.ka_params = self._blast.ka_params + (float(self._value),)
def read(handle, debug=0):
"""Returns a single Blast record (assumes just one query).
Uses the BlastParser internally.
This function is for use when there is one and only one BLAST
result in your XML file.
Use the Bio.Blast.NCBIXML.parse() function if you expect more than
one BLAST record (i.e. if you have more than one query sequence).
"""
iterator = parse(handle, debug)
try:
first = next(iterator)
except StopIteration:
first = None
if first is None:
raise ValueError("No records found in handle")
try:
second = next(iterator)
except StopIteration:
second = None
if second is not None:
raise ValueError("More than one record found in handle")
return first
def parse(handle, debug=0):
"""Returns an iterator a Blast record for each query.
Incremental parser, this is an iterator that returns
Blast records. It uses the BlastParser internally.
handle - file handle to and XML file to parse
debug - integer, amount of debug information to print
This is a generator function that returns multiple Blast records
objects - one for each query sequence given to blast. The file
is read incrementally, returning complete records as they are read
in.
Should cope with new BLAST 2.2.14+ which gives a single XML file
for multiple query records.
Should also cope with XML output from older versions BLAST which
gave multiple XML files concatenated together (giving a single file
which strictly speaking wasn't valid XML).
"""
from xml.parsers import expat
BLOCK = 1024
MARGIN = 10 # must be at least length of newline + XML start
XML_START = "<?xml"
text = handle.read(BLOCK)
pending = ""
if not text:
# NO DATA FOUND!
raise ValueError("Your XML file was empty")
while text:
# We are now starting a new XML file
if not text.startswith(XML_START):
raise ValueError("Your XML file did not start with %s... "
"but instead %s"
% (XML_START, repr(text[:20])))
expat_parser = expat.ParserCreate()
blast_parser = BlastParser(debug)
expat_parser.StartElementHandler = blast_parser.startElement
expat_parser.EndElementHandler = blast_parser.endElement
expat_parser.CharacterDataHandler = blast_parser.characters
expat_parser.Parse(text, False)
while blast_parser._records:
record = blast_parser._records[0]
blast_parser._records = blast_parser._records[1:]
yield record
while True:
# Read in another block of the file...
text, pending = pending + handle.read(BLOCK), ""
if not text:
# End of the file!
expat_parser.Parse("", True) # End of XML record
break
# Now read a little bit more so we can check for the
# start of another XML file...
pending = handle.read(MARGIN)
if ("\n" + XML_START) not in (text + pending):
# Good - still dealing with the same XML file
expat_parser.Parse(text, False)
while blast_parser._records:
yield blast_parser._records.pop(0)
else:
# This is output from pre 2.2.14 BLAST,
# one XML file for each query!
# Finish the old file:
text, pending = (text + pending).split("\n" + XML_START, 1)
pending = XML_START + pending
expat_parser.Parse(text, True) # End of XML record
while blast_parser._records:
yield blast_parser._records.pop(0)
# Now we are going to re-loop, reset the
# parsers and start reading the next XML file
text, pending = pending, ""
break
# this was added because it seems that the Jython expat parser
# was adding records later then the Python one
while blast_parser._records:
yield blast_parser._records.pop(0)
# At this point we have finished the first XML record.
# If the file is from an old version of blast, it may
# contain more XML records (check if text=="").
assert pending == ""
assert len(blast_parser._records) == 0
# We should have finished the file!
assert text == ""
assert pending == ""
assert len(blast_parser._records) == 0
if __name__ == '__main__':
import sys
with open(sys.argv[1]) as handle:
r_list = parse(handle)
for r in r_list:
# Small test
print('Blast of %s' % r.query)
print('Found %s alignments with a total of %s HSPs'
% (len(r.alignments),
reduce(lambda a, b: a + b,
[len(a.hsps) for a in r.alignments])))
for al in r.alignments:
print("%s %i bp %i HSPs" % (al.title[:50], al.length, len(al.hsps)))
# Cookbook example
E_VALUE_THRESH = 0.04
for alignment in r.alignments:
for hsp in alignment.hsps:
if hsp.expect < E_VALUE_THRESH:
print('*****')
print('sequence %s' % alignment.title)
print('length %i' % alignment.length)
print('e value %f' % hsp.expect)
print(hsp.query[:75] + '...')
print(hsp.match[:75] + '...')
print(hsp.sbjct[:75] + '...')