/Bio/Blast/NCBIStandalone.py
Python | 2199 lines | 2065 code | 37 blank | 97 comment | 54 complexity | c753c26957376830ac0068a5828c4d28 MD5 | raw file
Large files files are truncated, but you can click here to view the full file
- # 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 (OBSOLETE).
- 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
- prefered. 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).
- Classes:
- LowQualityBlastError Except that indicates low quality query sequences.
- BlastParser Parses output from blast.
- BlastErrorParser Parses output and tries to diagnose possible errors.
- PSIBlastParser Parses output from psi-blast.
- Iterator Iterates over a file of blast results.
- _Scanner Scans output from standalone BLAST.
- _BlastConsumer Consumes output from blast.
- _PSIBlastConsumer Consumes output from psi-blast.
- _HeaderConsumer Consumes header information.
- _DescriptionConsumer Consumes description information.
- _AlignmentConsumer Consumes alignment information.
- _HSPConsumer Consumes hsp information.
- _DatabaseReportConsumer Consumes database report information.
- _ParametersConsumer Consumes parameters information.
- Functions:
- blastall Execute blastall (OBSOLETE).
- blastpgp Execute blastpgp (OBSOLETE).
- rpsblast Execute rpsblast (OBSOLETE).
- For calling the BLAST command line tools, we encourage you to use the
- command line wrappers in Bio.Blast.Applications - the three functions
- blastall, blastpgp and rpsblast are considered to be obsolete now, and
- are likely to be deprecated and then removed in future releases.
- """
- import warnings
- warnings.warn("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.", PendingDeprecationWarning)
- import os
- import re
- import StringIO
- from Bio import File
- from Bio.ParserSupport import *
- from Bio.Blast import Record
- from Bio.Application import _escape_filename
- 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:
- """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 1:
- 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() : break
- if "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
- line.find('No hits found') == -1):
- 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 line.find("ERROR:") != -1 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, 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)
- 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 1:
- 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 1:
- 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 1:
- # 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, 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 1:
- 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, 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 line.find('Lambda') != -1:
- break
- read_and_call(uhandle, consumer.noevent, start='Lambda')
- read_and_call(uhandle, consumer.ka_params)
- #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, 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 ommitted 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:
- 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:
- 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:
- # 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,
- self._start_index+len(start))
- # subtract 1 for the space
- self._name_length = self._start_index - 1
- self._start_length = self._seq_index - self._start_index - 1
- self._seq_length = line.rfind(end) - self._seq_index - 1
-
- #self._seq_index = line.index(seq)
- ## subtract 1 for the space
- #self._seq_length = line.rfind(end) - self._seq_index - 1
- #self._start_index = line.index(start)
- #self._start_length = self._seq_index - self._start_index - 1
- #self._name_length = self._start_index
- # Extract the information from the line
- name = line[:self._name_length]
- name = name.rstrip()
- start = line[self._start_index:self._start_index+self._start_length]
- start = start.rstrip()
- if start:
- start = _safe_int(start)
- end = line[self._seq_index+self._seq_length:].rstrip()
- if end:
- end = _safe_int(end)
- seq = line[self._seq_index:self._seq_index+self._seq_length].rstrip()
- # right pad the sequence with spaces if necessary
- if len(seq) < self._seq_length:
- seq = seq + ' '*(self._seq_length-len(seq))
-
- # I need to make sure the sequence is aligned correctly with the query.
- # First, I will find the length of the query. Then, if necessary,
- # I will pad my current sequence with spaces so that they will line
- # up correctly.
- # Two possible things can happen:
- # QUERY
- # 504
- #
- # QUERY
- # 403
- #
- # Sequence 504 will need padding at the end. Since I won't know
- # this until the end of the alignment, this will be handled in
- # end_alignment.
- # Sequence 403 will need padding before being added to the alignment.
- align = self._multiple_alignment.alignment # for convenience
- align.append((name, start, seq, end))
- # This is old code that tried to line up all the sequences
- # in a multiple alignment by using the sequence title's as
- # identifiers. The problem with this is that BLAST assigns
- # different HSP's from the same sequence the same id. Thus,
- # in one alignment block, there may be multiple sequences with
- # the same id. I'm not sure how to handle this, so I'm not
- # going to.
-
- # # If the sequence is the query, then just add it.
- # if name == 'QUERY':
- # if len(align) == 0:
- # align.append((name, start, seq))
- # else:
- # aname, astart, aseq = align[0]
- # if name != aname:
- # raise ValueError, "Query is not the first sequence"
- # aseq = aseq + seq
- # align[0] = aname, astart, aseq
- # else:
- # if len(align) == 0:
- # raise ValueError, "I could not find the query sequence"
- # qname, qstart, qseq = align[0]
- #
- # # Now find my sequence in the multiple alignment.
- # for i in range(1, len(align)):
- # aname, astart, aseq = align[i]
- # if name == aname:
- # index = i
- # break
- # else:
- # # If I couldn't find it, then add a new one.
- # align.append((None, None, None))
- # index = len(align)-1
- # # Make sure to left-pad it.
- # aname, astart, aseq = name, start, ' '*(len(qseq)-len(seq))
- #
- # if len(qseq) != len(aseq) + len(seq):
- # # If my sequences are shorter than the query sequence,
- # # then I will need to pad some spaces to make them line up.
- # # Since I've already right padded seq, that means aseq
- # # must be too short.
- # aseq = aseq + ' '*(len(qseq)-len(aseq)-len(seq))
- # aseq = aseq + seq
- # if astart is None:
- # astart = start
- # align[index] = aname, astart, aseq
- def end_alignment(self):
- # Remove trailing newlines
- if self._alignment:
- self._alignment.title = self._alignment.title.rstrip()
- # This code is also obsolete. See note above.
- # If there's a multiple alignment, I will need to make sure
- # all the sequences are aligned. That is, I may need to
- # right-pad the sequences.
- # if self._multiple_alignment is not None:
- # align = self._multiple_alignment.alignment
- # seqlen = None
- # for i in range(len(align)):
- # name, start, seq = align[i]
- # if seqlen is None:
- # seqlen = len(seq)
- # else:
- # if len(seq) < seqlen:
- # seq = seq + ' '*(seqlen - len(seq))
- # align[i] = name, start, seq
- # elif len(seq) > seqlen:
- # raise ValueError, \
- # "Sequence %s is longer than the query" % name
-
- # Clean up some variables, if they exist.
- try:
- del self._seq_index
- del self._seq_length
- del self._start_index
- del self._start_length
- del self._name_length
- except AttributeError:
- pass
- class _HSPConsumer:
- def start_hsp(self):
- self._hsp = Record.HSP()
- def score(self, line):
- self._hsp.bits, self._hsp.score = _re_search(
- r"Score =\s*([0-9.e+]+) bits \(([0-9]+)\)", line,
- "I could not find the score in line\n%s" % line)
- self._hsp.score = _safe_float(self._hsp.score)
- self._hsp.bits = _safe_float(self._hsp.bits)
- x, y = _re_search(
- r"Expect\(?(\d*)\)? = +([0-9.e\-|\+]+)", line,
- "I could not find the expect in line\n%s" % line)
- if x:
- self._hsp.num_alignments = _safe_int(x)
- else:
- self._hsp.num_alignments = 1
- self._hsp.expect = _safe_float(y)
- def identities(self, line):
- x, y = _re_search(
- r"Identities = (\d+)\/(\d+)", line,
- "I could not find the identities in line\n%s" % line)
- self._hsp.identities = _safe_int(x), _safe_int(y)
- self._hsp.align_length = _safe_int(y)
- if line.find('Positives') != -1:
- x, y = _re_search(
- r"Positives = (\d+)\/(\d+)", line,
- "I could not find the positives in line\n%s" % line)
- self._hsp.positives …
Large files files are truncated, but you can click here to view the full file